/*************************************************************
File contains procedures:
  it_create_blocks - to allocate space for a block structure
  it_init_blocks - to initialize block structure
  it_assemble_blocks - to assemble entries to the global stiffness matrix
                           and the global load vector using the provided local 
                           stiffness matrix and load vector
  it_create_blocks_dia - to create preconditioner blocks corresponding
                       to small subdomains of neighboring elements
  it_factor_blocks_dia - to factorize the diagonal blocks of stiffness matrix
  it_compres - to compute the residual of the not preconditioned 
	system of equations, v = ( b - Ax )
  it_blsiter - to perform one iteration of block Gauss-Seidel
	or block Jacobi algorithm - for small blocks
  it_blliter - to perform one iteration of block Gauss-Seidel
	or block Jacobi algorithm - for large blocks
  it_rhsub - to perform forward reduction and back-substitution for ILU
           preconditioning

Only as pattern for future developments:
  it_mfaiter - to perform matrix vector multiplication (possibly in
	matrix free manner) and additive Schwarz approximate solve
  it_mfmiter - to perform one iteration of block Gauss-Seidel
 	for matrix-free GMRES

------------------------------  			
History:        
	02.2002 - Krzysztof Banas, initial version		
*************************************************************************/

#include<stdlib.h>
#include<stdio.h>
#include<string.h>
#include<math.h>

/* internal information for the solver module */
#include "./it_bliter.h"

/* LAPACK and BLAS procedures */
#include "lin_alg_intf.h"

/* minimal number of dofs to make use of DGEMV reasonable */
#define MIN_DOF_DGEMV 10000
#define SMALL 1e-15 /* approximation of round-off error */
#define TOL 1e-9    /* default tolerance for different iterative processes */

#define Max_dof_block_dia 1500
#define MAX_BL_BL 20

const int ITC_MAX_BL_NGB = 200; /* maximal no of block's neighbors */
const int ITC_MAX_BL_DIA_NGB = 200;/* maximal number of dia block's neighbors */


/*---------------------------------------------------------
it_create_blocks - to allocate space for a block structure
---------------------------------------------------------*/
int it_create_blocks(
  int Solver_id,   /* in: solver ID (used to identify the subproblem) */
  int Level_id,    /* in: level ID */
  int Nrblocks,    /* in: number of DOF blocks */
  int Nrdof_glob,  /* in: total number of DOFs */
  int* Nrdofbl,	   /* in: list of numbers of dofs in a block */
  int* Posglob,	   /* in: list of global numbers of first dof */
  int* Nroffbl,	   /* in: list of numbers of off diagonal blocks */
  // if Nroffbl[iblock]==0 - iblock is a ghost block with no rows associated with
  int** L_offbl	   /* in: list of lists of off diagonal blocks */
	)
{

/* local variables */
  itt_levels *it_level; /* mesh levels */
  itt_blocks *block;	/* to simplify typing */

/* auxiliary variable */
  int i, iblock, bloff, nrblocks, nrdofbl, nroffbl;

/*++++++++++++++++ executable statements ++++++++++++++++*/

  /* get pointer to the level structure */
  it_level = &itv_solver[Solver_id].level[Level_id];

  it_level->Nrblocks = Nrblocks;
  it_level->Nrdofgl = Nrdof_glob;

  /* allocate space for pointers to blocks */
  nrblocks = it_level->Nrblocks;
  it_level->Block = (itt_blocks **)malloc((nrblocks+1)*sizeof(itt_blocks *));
  if( it_level->Block==NULL ) {
    printf("Not enough memory for allocating block structure for stiffness matrix\n");
    exit(-1);
  }
    
  for(iblock=1;iblock<=nrblocks;iblock++){
      
    /* allocate space for a new block */
    it_level->Block[iblock] = (itt_blocks *)malloc(sizeof(itt_blocks));
    if(it_level->Block[iblock]==NULL){
      printf("it_create_blocks: no space for block %d\n", iblock);
      exit(1);
    }
      
    block = it_level->Block[iblock];
      
    /* specify block data */
    block->Ndof = Nrdofbl[iblock];
    block->Posg = Posglob[iblock];
    if(Nroffbl[iblock]>0){
      block->Lngb=it_ivector(Nroffbl[iblock]+1,"Lngb in it_create_blocks");
      block->Lngb[0] = Nroffbl[iblock];
      for(i=1;i<=Nroffbl[iblock];i++){
	block->Lngb[i] = L_offbl[iblock][i-1];
      }
    }
    else block->Lngb=NULL;
    // all ghost blocks are identified by block->Lngb==NULL
      
  }
    
  for(iblock=1;iblock<=nrblocks;iblock++){
      
    block = it_level->Block[iblock];
      
    nrdofbl = block->Ndof;


    /* for ordinary (not ghost) blocks */
    if(block->Lngb!=NULL){

      /* allocate space for diagonal matrix, rhs vector and */
      /* auxiliary ips vector for pivoting information */
      block->Dia=it_dvector(nrdofbl*nrdofbl,"Dia in it_create_blocks");
      block->Rhs=it_dvector(nrdofbl,"Rhs in it_create_blocks");
      block->Ips=it_ivector(nrdofbl,"ips it_create_blocks");
      
      /* array of pointers to off-diagonal blocks */
      block->Aux=(double **)malloc(block->Lngb[0]*sizeof(double *)); 
      
      for(i=1;i<=block->Lngb[0];i++) {
	
	bloff=block->Lngb[i];
	nroffbl = it_level->Block[bloff]->Ndof;
	
	/* allocate space for a new off-diagonal matrix */
	block->Aux[i-1]=it_dvector(nrdofbl*nroffbl,"Aux in it_create_blocks");
	
      }
    }
    else{

      block->Dia=NULL;
      block->Rhs=NULL;
      block->Ips=NULL;
      block->Aux=NULL;

    }
  }    

  return(0);
}

/*---------------------------------------------------------
  it_init_blocks - to initialize block structure
---------------------------------------------------------*/
int it_init_blocks(
  int Solver_id,   /* in: solver ID (used to identify the subproblem) */
  int Level_id,    /* in: level ID */
  int Comp_type    /* in: indicator for the scope of computations: */
                   /*   ITC_SOLVE - solve the system */
                   /*   ITC_RESOLVE - resolve for the new rhs vector */
  )

{

/* local variables */
  itt_levels *it_level; /* mesh levels */
  itt_blocks *block;	/* to simplify typing */

/* auxiliary variable */
  int i, iblock, ineig, bloff, nrblocks, nrdofbl, nroffbl;

/*++++++++++++++++ executable statements ++++++++++++++++*/

  it_level = &itv_solver[Solver_id].level[Level_id];

  /* allocate space for pointers to blocks */
  nrblocks = it_level->Nrblocks;
  for(iblock=1;iblock<=nrblocks;iblock++){

    block = it_level->Block[iblock];

    /* for ordinary (not ghost) blocks */
    if(block->Lngb!=NULL){
      
      nrdofbl = block->Ndof;
      
      /* initialize blocks forming RHS vector */
      it_d_zero(block->Rhs,nrdofbl);
      
      if(Comp_type==ITC_SOLVE){

	/* initialize arrays forming stiffness matrix */
	it_d_zero(block->Dia,nrdofbl*nrdofbl);

	for(ineig=1;ineig<=block->Lngb[0];ineig++) {
	
	  bloff=block->Lngb[ineig];
	  nroffbl = it_level->Block[bloff]->Ndof;
	
	  for(i=0;i<nrdofbl*nroffbl;i++) block->Aux[ineig-1][i]=0;
      
	}
      }
    }
  }

  return(0);
}


/*------------------------------------------------------------
  it_assemble_blocks - to assemble entries to the global stiffness matrix
                           and the global load vector using the provided local 
                           stiffness matrix and load vector
------------------------------------------------------------*/
int it_assemble_blocks( 
                         /* returns: >=0 - success code, <0 - error code */
  int Solver_id,         /* in: solver ID (used to identify the subproblem) */
  int Level_id,          /* in: level ID */
  int Comp_type,         /* in: indicator for the scope of computations: */
                         /*   ITC_SOLVE - solve the system */
                         /*   ITC_RESOLVE - resolve for the new rhs vector */
  int Nr_dof_bl,         /* in: number of global dof blocks */
                         /*     associated with the local stiffness matrix */
  int* L_bl_id,          /* in: list of dof blocks' IDs */
  int* L_bl_nrdof,       /* in: list of blocks' numbers of dof */
  double* Stiff_mat,     /* in: stiffness matrix stored columnwise */
  double* Rhs_vect,      /* in: rhs vector */
  char* Rewr_dofs         /* in: flag to rewrite or sum up entries */
                         /*   'T' - true, rewrite entries when assembling */
                         /*   'F' - false, sum up entries when assembling */
  )
{

  /* pointers to solver structure and stiffness matrix blocks' info */
  itt_levels *it_level; /* mesh levels */
  itt_blocks *block_i, *block_j;	/* to simplify typing */

  /* auxiliary variables */
  int iblock, jblock, nrdof_i, nrdof_j, posbl_i, posbl_j;
  int posloc_i, posloc_j, nrdof, nrdof_glob, ibl;
  int i,j,k, iaux, jaux;

/*kbw
  static FILE *fp=NULL;
/*kew*/

/*++++++++++++++++ executable statements ++++++++++++++++*/

/*kbw
  if(fp==NULL) fp=fopen("stiff_mat","a"); 
/*kew*/

  /* get pointer to the level structure */
  it_level = &itv_solver[Solver_id].level[Level_id];
  
  /* compute local stiffness matrix nrdof */
  nrdof=0;
  for(iblock=0;iblock<Nr_dof_bl;iblock++){
    nrdof += L_bl_nrdof[iblock];
  }
  
  /* loop over stiffness matrix blocks */
  posloc_i=0; 
  for(iblock=0;iblock<Nr_dof_bl;iblock++){

    block_i =  it_level->Block[L_bl_id[iblock]];

    /* number of dofs for a block */
    nrdof_i = L_bl_nrdof[iblock];

/*kbw
    if(Solver_id != -1 && block_i->Posg == 0){
      printf("\nin horizontal block_i %d, posglob_i %d, bl_id %d, nrdofs_i %d\n\n",
	     iblock, block_i->Posg, L_bl_id[iblock], L_bl_nrdof[iblock]);
      getchar();
    } 
/*kew*/

#ifdef DEBUG
    if(nrdof_i!=block_i->Ndof){
      printf("Error in block structure for SM in assemble!, i %d != %d\n",
	     nrdof_i, block_i->Ndof);
      exit(-1);
    }
#endif
            
    posloc_j=0;
    for(jblock=0;jblock<Nr_dof_bl;jblock++){
      
      block_j =  it_level->Block[L_bl_id[jblock]];
      
      /* number of dofs for a block */
      nrdof_j = L_bl_nrdof[jblock];
      
      
/*kbw
      if(Solver_id != -1 && block_i->Posg == 0){
	printf("block_i %d, posglob_i %d, nrdofs_i %d\n",
	       iblock, block_i->Posg, block_i->Ndof);
	printf("block_j %d, posglob_i %d, nrdofs_i %d\n",
	       jblock, block_j->Posg, block_j->Ndof);
      } 
/*kew*/

      /* for ordinary (not ghost) blocks */
      if(block_j->Lngb!=NULL){

#ifdef DEBUG
	if(nrdof_j!=block_j->Ndof){
	  printf("Error in block structure for SM in assemble!, j %d != %d\n",
		 nrdof_j, block_j->Ndof);
	  exit(-1);
	}
#endif
      
	iaux = 0;
	jaux = posloc_j + posloc_i*nrdof;
      
	if(block_i==block_j){
	
	  for(i=0;i<nrdof_i;i++){

/*kbw
	    if(Solver_id != -1 && block_i->Posg == 0){
	      printf("assembling Dia (%d) entries \nfrom local : %d - %d ", 
		     L_bl_id[iblock], jaux, jaux+nrdof_j-1 );
	      printf("to global %d - %d \nvalues: (local global): ", 
		     block_j->Posg + (block_i->Posg+i)*it_level->Nrdofgl,
		     block_j->Posg + (block_i->Posg+i)*it_level->Nrdofgl + nrdof_j-1);
	    } 
/*kew*/

	    for(j=0;j<nrdof_j;j++){
	      block_i->Dia[iaux+j] += Stiff_mat[jaux+j];	    
	    
/*kbw
	      if(Solver_id != -1 && block_i->Posg == 0){
		printf("%20.15lf   %20.15lf\n",
		       Stiff_mat[jaux+j], block_i->Dia[iaux+j]);
	      } 
/*kew*/
	  
	    }
	  
	    jaux += nrdof;
	    iaux += nrdof_j;
	  }
	
	  
	}
	else{
	
/*kbw 
if(Solver_id == 2 && block_i->Posg < 20){
	printf("block_dia %d, posglob %d, nrdof %d\n",
		L_bl_id[jblock], block_j->Posg, nrdof_j);
} 
/*kew*/
/*kbw if(Solver_id == 2 && block_i->Posg < 20){
	printf("block_off %d, posglob %d, nrdof %d\n",
		L_bl_id[iblock], block_i->Posg, nrdof_i);
} /*kew*/
	
	  ibl = it_chk_list(L_bl_id[iblock],&block_j->Lngb[1],block_j->Lngb[0]);
	
/*kbw if(Solver_id == 2 && block_i->Posg < 20){
	printf("assembling Aux[%d] - %d values:\n",
		ibl-1,nrdof_j*nrdof_i );
} /*kew*/
	
	  for(i=0;i<nrdof_i;i++){

/*kbw if(Solver_id == 2 && block_i->Posg < 20){
	printf("assembling Aux[%d] entries \nfrom local : %d - %d ", 
	       ibl-1, jaux, jaux+nrdof_j-1 );
	printf("to global %d - %d \nvalues: (local global): ", 
	       block_j->Posg + (block_i->Posg+i) * it_level->Nrdofgl,
	       block_j->Posg + (block_i->Posg+i)*it_level->Nrdofgl + nrdof_j-1 );
} /*kew*/

	    for(j=0;j<nrdof_j;j++){
	      block_j->Aux[ibl-1][iaux+j] += Stiff_mat[jaux+j];
	    
/*kbw if(Solver_id == 2 && block_i->Posg < 20){
	    printf("%20.15lf   %20.15lf\n",
		    Stiff_mat[jaux+j], block_j->Aux[ibl-1][iaux+j]);
} /*kew*/
	    
	    }
	  
	    jaux += nrdof;
	    iaux += nrdof_j;
	  }
	
	}

	/*kbw if(Solver_id == 2 && block_i->Posg < 20){
	printf("\n");
	} /*kew*/

      } /* end if ordinary (not ghost) jblock */
      
      posloc_j += nrdof_j;
    }
    
    /* for assembling RHSV loop over iblocks should be treated as loop over */
    /* VERTICAL blocks, i.e. now iblock is a subsequent vertical block */

    /* for ordinary (not ghost) blocks */
    if(block_i->Lngb!=NULL){
      
      /*kbw if(Solver_id == 2 && block_i->Posg < 20){
      printf("assembling Rhs (%d) entries \nfrom local : %d - %d ",
	     L_bl_id[iblock], posloc_i, posloc_i+nrdof_i-1 );
      printf("to global %d - %d \nvalues: (local global): ", 
	     block_i->Posg, block_i->Posg+nrdof_i-1 );
      } /*kew*/
      
      for(i=0;i<nrdof_i;i++) {
	
	/* assemble right hand side block's entries */
	block_i->Rhs[i] += Rhs_vect[posloc_i+i];
	
	/*kbw if(Solver_id == 2 && block_i->Posg < 20){
	printf("%20.15lf   %20.15lf\n",
	       Rhs_vect[posloc_i+i], block_i->Rhs[i]);
	} /*kew*/
	
      }
      /*kbw if(Solver_id == 2 && block_i->Posg < 20){
      printf("\n");
      } /*kew*/

    }

      /*kbw if(Solver_id == 2 && block_i->Posg < 20){
  getchar();
      } /*kew*/

    posloc_i += nrdof_i;
  }
  

  return(1);
}

/*---------------------------------------------------------
it_create_blocks_dia - to create preconditioner blocks corresponding
                       to small subdomains of neighboring elements
---------------------------------------------------------*/
int it_create_blocks_dia( /* returns:   >0 number of diagonal blocks */
                          /*	       <=0 - error */
  int Solver_id,   /* in: solver ID (used to identify the subproblem) */
  int Level_id    /* in: level ID */
  )
{

  /* pointers to solver structure and stiffness matrix blocks' info */
  itt_levels *it_level; /* mesh levels */

/* local variables */
  itt_blocks_dia *block_dia;/* to simplify typing */
  int iblock, ib_dia;	/* counters for blocks */
  int nrdof_dia;		/* number of dofs in a block */
  int nmfa;	        /* number of faces */
  
  int* tags;		/* to tag objects for which blocks were created */
  int nr_neighb, bl_neighb[100];/* number and list of block's neighbors */
  int **temp_list, i_temp_list;	/* temporary list for block_dia blocks */
  int nr_bl_block;	/* number of elementary blocks in a dia block */
  int neig_loc[ITC_MAX_BL_DIA_NGB]; /* temp array with dia block neighbors */
  int **big_block; /* temporary array with big block numbers for small blocks */


/* auxiliary variables */
  int i, j, k, l, iaux, jaux, kaux, laux, ibl, nrnos, nrneig, ineig;

/*++++++++++++++++ executable statements ++++++++++++++++*/

  /* get pointer to the level structure */
  it_level = &itv_solver[Solver_id].level[Level_id];

  if(abs(it_level->Block_type)==1){

/* allocate space for pointers to blocks */
    it_level->Block_dia = 
      (itt_blocks_dia **)malloc((it_level->Nrblocks+1)
			       *sizeof(itt_blocks_dia *)); 
  
/* dia->Dia blocks are just inverted Dia blocks */
    it_level->Nrblocks_dia = 0;
    ib_dia = 0;
    for(iblock=1;iblock<=it_level->Nrblocks;iblock++){

      /* for ordinary (not ghost) blocks */
      if(it_level->Block[iblock]->Lngb!=NULL){

	it_level->Nrblocks_dia++;
	ib_dia++;

	/* allocate space for block_dia structure */
	it_level->Block_dia[ib_dia] = 
	  (itt_blocks_dia *)malloc(sizeof(itt_blocks_dia));
	if(it_level->Block_dia[ib_dia]==NULL){
	  printf("it_create_blocks_dia: no space for block %d\n", ib_dia);
	  return(1);
	}

	/* number of small blocks (nodes) in dia is equal to 1*/
	nrnos=1;
    
	/* prepare diagonal block for assembling */
	block_dia = it_level->Block_dia[ib_dia];
	block_dia->Lsmall=it_ivector(nrnos+1,"block_dia->Lsmall in create_blocks_dia");

	if(it_level->Solver==GMRES&&it_level->GMRES_type==MATRIX_FREE)
	  block_dia->Lelem=it_ivector(ITC_MAX_BL_DIA_NGB+1,
				      "block_dia->Lelem in create_blocks_dia");
	else block_dia->Lelem=NULL;

	block_dia->Lpos=it_ivector(nrnos+1,
				   "block_dia->Lpos in create_blocks_dia");
      
	/*initialize*/
	if(it_level->Solver==GMRES&&it_level->GMRES_type==MATRIX_FREE)
	  for(ibl=0;ibl<=ITC_MAX_BL_DIA_NGB;ibl++) block_dia->Lelem[ibl]=0;

	/* number of small blocks */
	block_dia->Lsmall[0] = nrnos;

	/* put small block number */
	block_dia->Lsmall[1] = iblock;

	/* number of neighbors and their numbers */
	if( it_level->Precon==MULTI_ILU||it_level->Precon==BLOCK_ILU ){
	  iaux = it_level->Block[iblock]->Lngb[0];
	  nrneig=0;
	  for(i=1;i<=iaux;i++){
	    jaux = it_level->Block[iblock]->Lngb[i];
	    if(it_level->Block[jaux]->Lngb!=NULL) nrneig++;
	  }
	  block_dia->Lneig=it_ivector(nrneig+1,"block_dia->Lneig in create_blocks_dia");
	  block_dia->Lneig[0] = nrneig;
	  nrneig=0;
	  for(i=1;i<=iaux;i++){
	    jaux = it_level->Block[iblock]->Lngb[i];
	    if(it_level->Block[jaux]->Lngb!=NULL) {
	      nrneig++;
	      block_dia->Lneig[nrneig]=it_level->Block[iblock]->Lngb[i];
	    }
	  }
	}
	else block_dia->Lneig=NULL;

	/* indicate local position of first dof in dia subarray */
	block_dia->Lpos[1] = 0;

	/* total number of dofs for dia */
	nrdof_dia = it_level->Block[iblock]->Ndof;

	/* total number of dofs for dia */
	block_dia->Lpos[0] = nrdof_dia;

	/* allocate space for dia and Ips */
	block_dia->Dia=
	  it_dvector(nrdof_dia*nrdof_dia,"Dia in create_blocks_dia");
	block_dia->Ips=
	  it_ivector(nrdof_dia,"Ips in create_blocks_dia");

      } /* if ordinary (not ghost) block */

    } /* for all diagonal blocks */

    /* for matrix free version create a list of block_dia integration entities */
    if(it_level->Solver==GMRES&&it_level->GMRES_type==MATRIX_FREE){
      printf("TBD - for matrix free create a list of block_dia blocks\n");
      printf(" max size < ITC_MAX_BL_DIA_NGB ???");
    }
    
  }
  else if(abs(it_level->Block_type)>1){ /* for larger blocks */

/* subdomains composed of blocks and its neighbors */
    if(    abs(it_level->Block_type)==2 || abs(it_level->Block_type)==3
	|| abs(it_level->Block_type)==4 || abs(it_level->Block_type)==5 
	|| abs(it_level->Block_type)==6 || abs(it_level->Block_type)==7 ){

      temp_list = it_imatrix(it_level->Nrblocks+1,ITC_MAX_BL_NGB+1,
			     "temp_list in create_block_dia");

/* for neighborhood blocks with different overlap */

/* tagged block are excluded from being internal block_dia blocks */
      tags = it_ivector(it_level->Nrblocks+1,"tags in create_blocks_dia");
      for(ibl=1;ibl<=it_level->Nrblocks;ibl++) tags[ibl] = 0;

/* allocate space for diagonal blocks */
      it_level->Nrblocks_dia = 0;

/* in a loop over all blocks */
      for(ibl=1;ibl<=it_level->Nrblocks;ibl++){

/* if block ordinary (not ghost) and not tagged */
	if( it_level->Block[ibl]->Lngb != NULL && (tags[ibl]==0 ||
			 abs(it_level->Block_type)==2 ||
			 abs(it_level->Block_type)==3 )){

/* increase the number of diagonal blocks and tag a block */
	  it_level->Nrblocks_dia++;
	  tags[ibl]=it_level->Nrblocks_dia;

/* initialize the counter for block_dia blocks */
	  nr_bl_block = 0;

/* put block on temporary list of block_dia blocks  */ 
	  nr_bl_block++;
	  temp_list[it_level->Nrblocks_dia][nr_bl_block]=ibl;

/* loop over block's neighbors */
/* i.e. other blocks connected by entries in the global stiffness matrix */
	  for(iaux=0;iaux<it_level->Block[ibl]->Lngb[0];iaux++){

	    ineig = it_level->Block[ibl]->Lngb[iaux+1];

/*kbw
printf("in neig %d (%d), LNGB %u, tag %d\n",
       ineig,iaux,it_level->Block[ineig]->Lngb,tags[ineig]);
/*kew*/

/* if ordinary (not ghost) block and (not tagged yet or we include overlap) */
	    if(  it_level->Block[ineig]->Lngb != NULL && (tags[ineig]==0 ||
			 abs(it_level->Block_type)==4 ||
			 abs(it_level->Block_type)==5 ) ){
      
/* put on temporary list of block_dia blocks  */ 
	      nr_bl_block++;
	      temp_list[it_level->Nrblocks_dia][nr_bl_block]=ineig;

/*kbw
printf("in neig %d (%d), put as block %d in block %d\n",
	ineig,iaux,nr_bl_block,it_level->Nrblocks_dia);
/*kew*/

/* tag */
	      if(abs(it_level->Block_type)==3 || 
		 abs(it_level->Block_type)==5 || 
		 abs(it_level->Block_type)==7 )
		tags[ineig] = it_level->Nrblocks_dia; 
	      
	      
	    } /* if neighbor included */
	    
	  }  /* for all block's neighbors: iaux */

/* the total number of blocks in a diagonal block */
	  if(nr_bl_block>=ITC_MAX_BL_NGB){
	    printf("Too much blocks in blockary block in it_create_blocks_dia\n");
	    printf("Change parameters for block creation or increase ITC_MAX_BL_NGB\n");
	    return(-1);
	  }
	  temp_list[it_level->Nrblocks_dia][0] = nr_bl_block;
	  
/*kbw
printf("Created block %d: %d blocks\n",
	it_level->Nrblocks_dia,nr_bl_block);
for(iaux=1;iaux<=nr_bl_block;iaux++) 
	printf("%d   ",temp_list[it_level->Nrblocks_dia][iaux]);
printf("\n");
 getchar();
/*kew*/
	  
	  
	} /* if block ordinary and not tagged  */
	
      } /* for all active blocks: ibl */
      
/*kbw
printf("Total number of blocks %d\n",it_level->Nrblocks_dia);
kew*/
 
     

/* allocate space for pointers to blocks */
      it_level->Block_dia = 
	(itt_blocks_dia **)malloc((it_level->Nrblocks_dia+1)
				 *sizeof(itt_blocks_dia *)); 
      
/* loop over diagonal blocks */
      for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){

/* allocate space for block_dia structure */
	it_level->Block_dia[ib_dia] = 
	  (itt_blocks_dia *)malloc(sizeof(itt_blocks_dia));
	if(it_level->Block_dia[ib_dia]==NULL){
          printf("it_level gmres: no space for block %d\n", ib_dia);
          exit(1);
	}

/* prepare diagonal block for assembling */
	block_dia = it_level->Block_dia[ib_dia];

/* renumbering: i_temp_list = it_level->Nrblocks_dia + 1 - ib_dia;  */
	i_temp_list = ib_dia;
      
 /* number of small blocks (blocks) in dia */
	iaux = temp_list[i_temp_list][0]; 
	block_dia->Lsmall=it_ivector(iaux+1,"block_dia in create_blocks_dia");
	block_dia->Lpos=it_ivector(iaux+1,"block_dia in create_blocks_dia");
	block_dia->Lelem=NULL;
	
/* number of small blocks */
	block_dia->Lsmall[0] = iaux;
	
/* compute the number of dofs for diagonal block */
	nrdof_dia = 0;

/* for each block in block_dia */
	for(ibl=1;ibl<=temp_list[i_temp_list][0];ibl++){

/* get small block number */
	  iblock=temp_list[i_temp_list][ibl];

/* put small block number */
	  block_dia->Lsmall[ibl] = iblock;

/* indicate local position of first dof in dia subarray */
	  block_dia->Lpos[ibl] = nrdof_dia;

/* update total number of dofs for dia */
	  nrdof_dia += it_level->Block[iblock]->Ndof;

/*kbw
printf("counting dofs of el %d (loc %d)\n",
	temp_list[i_temp_list][ibl],ibl);
printf("small block %d, nrdof %d, total nrdof %d\n",
	iblock,it_level->Block[iblock]->Ndof,nrdof_dia);
/*kew*/


	} /* for all  blocks in a block */


/* total number of dofs for dia */
	block_dia->Lpos[0] = nrdof_dia;

/* allocate space for dia and Ips */
	block_dia->Dia=
	  it_dvector(nrdof_dia*nrdof_dia,"Dia in create_blocks_dia");
	block_dia->Ips=
	  it_ivector(nrdof_dia,"Ips in create_blocks_dia");


/*kbw
  printf("In block %d, %d small blocks (blocks):\n",
	 ib_dia,iaux);
 for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){
   printf("%d (%d)  ",block_dia->Lsmall[iblock],
	  it_level->Mesh_elems[block_dia->Lsmall[iblock]]);
 }
 printf("\n");
/*kew*/

/*kbw
printf("allocated %d dofs in block %d\n",
	nrdof_dia*nrdof_dia, ib_dia);
getchar();
kew*/


      } /* for all diagonal blocks: ib_dia */
      
      for(i=0;i<=it_level->Nrblocks;i++) free(temp_list[i]);
      free(temp_list);
      free(tags);
      
    } /* if specific block type */
   
    
/* for block ILU factorization */
    if( it_level->Precon==MULTI_ILU ||  it_level->Precon==BLOCK_ILU ){
      big_block = it_imatrix(it_level->Nrblocks+1,MAX_BL_BL,"in create_dia");
      for(iblock=1;iblock<=it_level->Nrblocks;iblock++) {
	for(i=0;i<MAX_BL_BL;i++) big_block[iblock][i]=0;
      }
      for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){
	block_dia=it_level->Block_dia[ib_dia];
	for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){
	  j=block_dia->Lsmall[iblock];
	  i=it_put_list(ib_dia, &big_block[j][1], MAX_BL_BL);
	  if(!i) {
	    big_block[j][0]++; 
	    if(big_block[j][0]>=MAX_BL_BL){
	      printf("Small block %d in more than MAX_BL_BL dia_blocks!",j);
	      printf("Change parameter MAX_BL_BL in it_bliter_blocks.c\n");
	      return(-1);
	    }
	  }

	}
      }
/*kbw
      for(iblock=1;iblock<=it_level->Nrblocks;iblock++) {
	printf("Small block %d, Big blocks:\n",iblock);
	for(i=1;i<=big_block[iblock][0];i++) printf("%d ",big_block[iblock][i]);
	printf("\n");
      }	
/*kew*/

      for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){
	block_dia=it_level->Block_dia[ib_dia];
	nrneig=0;
	for(i=0;i<ITC_MAX_BL_DIA_NGB;i++) neig_loc[i]=0;
	for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){
	  iaux=block_dia->Lsmall[iblock];	  
	  for(i=1;i<=big_block[iaux][0];i++){
	    k=big_block[iaux][i];

/*kbw
	    printf(" block %d, small block %d , big_block %d \n",
		   ib_dia, iaux, k);
/*kew*/

	    if(k!=ib_dia){
	      l=it_put_list(k, neig_loc, ITC_MAX_BL_DIA_NGB );
	      if(!l) {
		nrneig++;
/*kbw
		printf("Block %d added as %d neighbor\n",k,nrneig);
		printf("Block %d added as %d neighbor\nneighbors:",k,nrneig);
		for(laux=0;laux<nrneig;laux++) {
		  printf("%d ",neig_loc[laux]);
		}
		printf("\n");
/*kew*/
	      }
	      if(l<0){
		printf("Too much neighbors (%d) for big block %d\n", 
		       nrneig, ib_dia);
		printf("Increase ITC_MAX_BL_DIA_NGB in create_blocks_dia\n");
		exit(1);
	      }
	    }
	  }
	  for(j=1;j<=it_level->Block[iaux]->Lngb[0];j++){
	    jaux=it_level->Block[iaux]->Lngb[j];
	    /* for parallel version only local blocks can be considered */
	    if(it_level->Block[jaux]->Lngb!=NULL){
	      for(i=1;i<=big_block[jaux][0];i++){
		k=big_block[jaux][i];

/*kbw
	      printf(" block %d, small block %d, big_block %d\n",
		     ib_dia, jaux, k);
/*kew*/

		if(k!=ib_dia){
		  l=it_put_list(k, neig_loc, ITC_MAX_BL_DIA_NGB );
		  if(!l) {
		    nrneig++;
/*kbw
		  printf("Block %d added as %d neighbor\nneighbors:",k,nrneig);
		  for(laux=0;laux<nrneig;laux++) {
		    printf("%d ",neig_loc[laux]);
		  }
		  printf("\n");
/*kew*/
		  }
		  if(l<0){
		    printf("Too much neighbors (%d) for big block %d\n", 
			   nrneig, ib_dia);
		    printf("Increase ITC_MAX_BL_DIA_NGB in create_blocks_dia\n");
		    exit(1);
		  }
		}
	      }
	    }
	  }
	}
	block_dia->Lneig=it_ivector(nrneig+1,
				    "block_dia->Lneig in create_blocks_dia");
	block_dia->Lneig[0]=nrneig;
	for(i=0;i<nrneig;i++) {
	  block_dia->Lneig[i+1]=neig_loc[i];
	}
      }

      for(i=0;i<=it_level->Nrblocks;i++) free(big_block[i]);
      free(big_block);
      
    } /* if block ILU preconditioner */
    else {
      for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){
	block_dia=it_level->Block_dia[ib_dia];
	block_dia->Lneig=NULL;
      }
    }

  } /* if larger blocks */
  
/*kbw
  for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){
    block_dia=it_level->Block_dia[ib_dia];
    printf("Big block %d, Neighbors:",ib_dia);
    for(i=1;i<=block_dia->Lneig[0];i++) printf("%d ",block_dia->Lneig[i]);
    printf("\n");
  }	
/*kew*/

/* allocate space for arrays of doubles: Aux[i] */ 
  if( it_level->Precon==MULTI_ILU||it_level->Precon==BLOCK_ILU ){
    for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){

      block_dia = it_level->Block_dia[ib_dia];
      nrdof_dia = block_dia->Lpos[0];

      nrneig = block_dia->Lneig[0];
      block_dia->Aux = (double **)malloc(nrneig*sizeof(double *)); 

      for(i=0;i<nrneig;i++){
	j = block_dia->Lneig[i+1];

	kaux = it_level->Block_dia[j]->Lpos[0];

/*kbw
  if(ib_dia>0){
    printf("Big block %d, nr_dof %d neighbor %d, nr_dofs %d\n",
	   ib_dia,nrdof_dia,j,kaux);
  getchar();
  }
/*kew*/

	block_dia->Aux[i] =
	  it_dvector(nrdof_dia*kaux,"Aux[] in create_blocks_dia");

      }

    } /* for all diagonal blocks */
  } /* if ILU preconditioner */ 

/*kbw
  printf("Order of blocks:\n");
  for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){
    block_dia=it_level->Block_dia[ib_dia];
    for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){
      iaux=block_dia->Lsmall[iblock];	
      printf("%d ",it_level->Mesh_elems[iaux]);
    }
 }
 printf("\n");
 getchar();
/*kew*/

  return(it_level->Nrblocks_dia);
}

/*---------------------------------------------------------
  it_factor_blocks_dia - to factorize the stiffness matrix, either only
                diagonal blocks or block ILU(0)
---------------------------------------------------------*/
int it_factor_blocks_dia( 
  int Solver_id,         /* in: solver ID (used to identify the subproblem) */
  int Level_id          /* in: level ID */
	)
{

  /* pointers to solver structure and stiffness matrix blocks' info */
  itt_levels *it_level; /* mesh levels */

/* local variables */
  int ilev;             /* current level number */
  itt_blocks *block;  /* to simplify typing */
  itt_blocks_dia *block_dia, *dia_neig, *block_k;/* to simplify typing */
  int iblock, ib_dia, kb_dia, jb_dia, ibs;	/* counters for blocks */
  int nrdofbl, nrdof_dia;	/* number of dofs in a block */
  int nrngbl;	    /* number of off diagonal small blocks (small neighbors) */
  int nrneig;	    /* number of off diagonal big blocks (big neighbors) */
  int ndof1,ndof2;	/* numbers of dofs for small blocks */
  int posloc1,posloc2;	/* positions in big DIA for small blocks */
  int ioff,koff;	/* off-diagonal block numbers for assembling DIA */

/* constants */
  int ione = 1;
  double done = 1.0;
  double dmone = -1.0;
  int lwork = 100000;

/* auxiliary variables */
  int i,j,k,iaux,jaux,kaux,ineig,jneig,ismneig;
  double daux,eaux,faux;
  //double *a_work;

/*++++++++++++++++ executable statements ++++++++++++++++*/

  /* get pointer to the level structure */
  it_level = &itv_solver[Solver_id].level[Level_id];

  /*kbw
  printf("before factorizing blocks_dia\n");
  for(iblock=1;iblock<=it_level->Nrblocks;iblock++){
    jaux=it_level->Block[iblock]->Ndof;
    printf("In block %d (%d rows, the first %d), %d blocks\n", iblock, 
	   it_level->Block[iblock]->Ndof, it_level->Block[iblock]->Posg,
	   it_level->Block[iblock]->Lngb[0]+1);
    if(it_level->Block[iblock]->Lngb!=NULL){
      
      printf("diagonal block (displayed by rows)\n");
      for(j=0;j<jaux;j++){
	for(k=0;k<jaux;k++){
	  printf("%12.3le", it_level->Block[iblock]->Dia[j+jaux*k]);
	}
	printf("\n");
      }
      
      iaux = it_level->Block[iblock]->Lngb[0];
      for(i=1;i<=iaux;i++){
	kaux = it_level->Block[iblock]->Lngb[i];
	printf("for block %d (from column %d) (entries displayed by rows!)\n",
	       kaux, it_level->Block[kaux]->Posg);
	for(j=0;j<jaux;j++){
	  for(k=0;k<jaux;k++){
	    printf("%12.3le", it_level->Block[iblock]->Aux[i-1][j+jaux*k]);
	  }
	printf("\n");
	}
      }
      
      printf("RHS\n");
      for(j=0;j<jaux;j++){
	printf("%12.3le", it_level->Block[iblock]->Rhs[j]);
      }
	printf("\n");
    }
    getchar();
  }
  /*kew*/

/*according to the value of the global variable it_level->Block_type */
  if(abs(it_level->Block_type)==0){
/* prepare block arrays */
    for(iblock=1;iblock<=it_level->Nrblocks;iblock++){
  
      if(it_level->Block[iblock]!=NULL&&it_level->Block[iblock]->Lngb!=NULL) {

        block = it_level->Block[iblock];
        nrdofbl = block->Ndof;
        nrngbl = block->Lngb[0];

/* for standard iterations as preconditioner prepare inverted arrays */
        if(nrdofbl==1){
          block->Rhs[0]=block->Rhs[0]/block->Dia[0];
          for(i=0;i<nrngbl;i++)block->Aux[i][0]=block->Aux[i][0]/block->Dia[0];
        }
        else{
/* LU decompose block->Dia */

          dgetrf_(&nrdofbl,&nrdofbl,block->Dia,&nrdofbl,block->Ips,&iaux);

/* resolve to find  block->Dia^(-1) * block->Rhs */
          kaux=1;

          dgetrs_("N",&nrdofbl,&kaux,block->Dia,&nrdofbl,block->Ips,
		block->Rhs,&nrdofbl,&iaux);

          for(i=0;i<nrngbl;i++){

/* find the number of dofs for the neighbor */
            j = block->Lngb[i+1];
            kaux = it_level->Block[j]->Ndof;

/* resolve to find  block->Dia^(-1) * block->Aux[i] */
            dgetrs_("N",&nrdofbl,&kaux,block->Dia,&nrdofbl,block->Ips,
		block->Aux[i],&nrdofbl,&iaux);

          }
        }

      } /* if block active */
    } /* for each block: iblock */

  } /* if Block_type=size == 0 (iterations with no dia blocks) */
  else { /* all other kinds of blocks */

    if(it_level->Precon!=MULTI_ILU&&it_level->Precon!=BLOCK_ILU){

/* for each diagonal block */
      for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){

/* if block is active */
	if(it_level->Block_dia[ib_dia]!=NULL){

/* to simplify typing */
	  block_dia=it_level->Block_dia[ib_dia];

/* assemble diagonal and suitable off diagonal submatrices from
   small blocks */
/* total number of dofs for DIA */
	  nrdof_dia = block_dia->Lpos[0];

/* initialize Dia and ips matrices */
	  for(i=0;i<nrdof_dia*nrdof_dia;i++) block_dia->Dia[i]=0.0;
	  for(i=0;i<nrdof_dia;i++) block_dia->Ips[i]=0;

/* loop over small blocks */
	  for(ibs=1;ibs<=block_dia->Lsmall[0];ibs++){

/* get into small block */
	    block=it_level->Block[block_dia->Lsmall[ibs]];

/* number of dofs for small block */
	    ndof1 = block->Ndof;

/* position of first dof in diagonal submatrix */
	    posloc1 = block_dia->Lpos[ibs];

/* assemble diagonal subarray */
	    for(i=0;i<ndof1;i++){
	      for(j=0;j<ndof1;j++){
		block_dia->Dia[posloc1+i+nrdof_dia*(posloc1+j)]
		  = block->Dia[i+ndof1*j];
	      }
	    }

/* if more than one small block */
	    if(block_dia->Lsmall[0]>1){

/* loop over small blocks */
	      for(iaux=1;iaux<=block_dia->Lsmall[0];iaux++){

/* if different from block ibs  */
		if(iaux!=ibs){

		  ioff = it_chk_list(block_dia->Lsmall[iaux],
				     &block->Lngb[1],block->Lngb[0]);

/* if there is an off-diagonal block */
		  if(ioff){

/* number of dofs for the off-diagonal submatrix */
		    ndof2 = it_level->Block[block_dia->Lsmall[iaux]]->Ndof;

/* position for the off-diagonal submatrix in dia */
		    posloc2 = block_dia->Lpos[iaux];

/* assemble off-diagonal subarray */
		    for(i=0;i<ndof1;i++){
		      for(j=0;j<ndof2;j++){
			block_dia->Dia[posloc1+i+nrdof_dia*(posloc2+j)]
			  = block->Aux[ioff-1][i+ndof1*j];
		      }
		    }

		  } /* if there is off-diagonal block */ 
		  
		} /* if two different blocks: iaux != ibs */

	      } /* for all small blocks: iaux */

	    } /* if more small blocks */ 

        } /* for all small blocks: ibs */

/* decompose Dia */

	  if(nrdof_dia==1) block_dia->Dia[0]= 1/block_dia->Dia[0];
	  else{

/* LU decompose block->Dia */
	    dgetrf_(&nrdof_dia,&nrdof_dia,block_dia->Dia,
		    &nrdof_dia,block_dia->Ips,&iaux);

	  }

	} /* if block is active */

      } /* for all blocks */

    } /* end if standard iterations preconditioning */
    else if(it_level->Precon==MULTI_ILU||it_level->Precon==BLOCK_ILU){

      /* prepare data for forward reduction and back substitution */
      printf("Preparing data for iterations\n");
  
/*       /\* loop over all big diagonal blocks *\/ */
/* #pragma omp parallel for  private(block_dia,jb_dia,ioff) */
/*       for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){ */
	
/* 	/\* to simplify typing *\/ */
/* 	block_dia=it_level->Block_dia[ib_dia]; */
	
/* 	block_dia->Llowerneig = (int *) malloc (200*sizeof(int)); */
/* 	block_dia->Llowerneig[0] = 0; */


/* 	//for all neighbours */
/* 	int iaux_neig; */
/* 	for(iaux_neig=1;iaux_neig<=block_dia->Lneig[0];iaux_neig++){  */

/* 	  jb_dia = block_dia->Lneig[iaux_neig]; */

/* 	  ioff=iaux_neig; */

/* 	  if(jb_dia<ib_dia){ */

/* 	/\* /\\* for all factorized (lower subdiagonal) entries of A *\\/ *\/ */
/* 	/\* for(jb_dia=1;jb_dia<ib_dia;jb_dia++){ *\/ */
	  
/* 	/\*   /\\* if two blocks are neighbors *\\/ *\/ */
/* 	/\*   ioff = it_chk_list(jb_dia,&block_dia->Lneig[1],block_dia->Lneig[0]); *\/ */
/* 	/\*   if(ioff){ *\/ */
	    
/* #ifdef DEBUG */
/* 	    if(block_dia->Lneig[ioff]!=jb_dia){ */
/* 	  printf("error in reading list of lower neighbors in decompose!\n"); */
/* 	  printf("ib_dia %d, jneig %d, jb_dia %d, ioff %d, ? %d\n", */
/*       ib_dia, block_dia->Llowerneig[0]+1, jb_dia, ioff, block_dia->Lneig[ioff]); */
/* 	      exit(1); */
/* 	    } */
/* 	    if(block_dia->Llowerneig[0]>=99){ */
/* 	      printf("increase the number of lower subdiagonal neighbors!\n"); */
/* 	      exit(1); */
/* 	    } */
/* #endif */
	    
/* 	    block_dia->Llowerneig[0]++; */
/* 	    //block_dia->Llowerneig[2*block_dia->Llowerneig[0]-1] = jb_dia; */
/* 	    //block_dia->Llowerneig[2*block_dia->Llowerneig[0]] = ioff; */

/* 	    // we must sort an array */
/* 	    i=2*block_dia->Llowerneig[0]-1; */
/* 	    while(i>0 && jb_dia<block_dia->Llowerneig[i-2]){ */
/* 	      block_dia->Llowerneig[i] = block_dia->Llowerneig[i-2]; */
/* 	      block_dia->Llowerneig[i+1] = block_dia->Llowerneig[i-1]; */
/* 	      i -= 2; */
/* 	    } */
/* 	    block_dia->Llowerneig[i] = jb_dia; */
/* 	    block_dia->Llowerneig[i+1] = ioff; */
	    
	    
/* 	  } */
/* 	} */
/*       } */
            
      /* loop over all big diagonal blocks */
#pragma omp parallel for default(none) private(block_dia,jb_dia,ioff,i) shared(it_level)
      //for(ib_dia=it_level->Nrblocks_dia;ib_dia>=1;ib_dia--){
      for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){

/* if block not active */
	if(it_level->Block_dia[ib_dia]==NULL){
	  printf("Block %d not active in it_factor_block_dia\n",ib_dia);
	  exit(0);
	}

	
	/* to simplify typing */
	block_dia=it_level->Block_dia[ib_dia];
	
	block_dia->Lupperneig = (int *) malloc (200*sizeof(int));
	block_dia->Lupperneig[0] = 0;

	block_dia->Llowerneig = (int *) malloc (200*sizeof(int));
	block_dia->Llowerneig[0] = 0;

	//for all neighbours
	int iaux_neig;
	for(iaux_neig=1;iaux_neig<=block_dia->Lneig[0];iaux_neig++){

	  jb_dia = block_dia->Lneig[iaux_neig];

	  ioff=iaux_neig;

	  if(jb_dia>ib_dia){


	/* for all upper diagonal blocks */
	//for(jb_dia=ib_dia+1;jb_dia<=it_level->Nrblocks_dia;jb_dia++){
	  
	  /* if two blocks are neighbors */
	  //ioff = it_chk_list(jb_dia,&block_dia->Lneig[1],block_dia->Lneig[0]);

	  //if(ioff){
	    
#ifdef DEBUG
	    if(block_dia->Lneig[ioff]!=jb_dia){
	  printf("error in reading list of upper neighbors in decompose!\n");
	  printf("ib_dia %d, jneig %d, jb_dia %d, ioff %d, ? %d\n",
      ib_dia, block_dia->Lupperneig[0]+1, jb_dia, ioff, block_dia->Lneig[ioff]);
	      exit(1);
	    }
	    if(block_dia->Lupperneig[0]>=99){
	      printf("increase the number of upper subdiagonal neighbors!\n");
	      exit(1);
	    }
#endif
	    
	    block_dia->Lupperneig[0]++;
	    //block_dia->Lupperneig[2*block_dia->Lupperneig[0]-1] = jb_dia;
	    //block_dia->Lupperneig[2*block_dia->Lupperneig[0]] = ioff;

	    // we must sort an array
	    i=2*block_dia->Lupperneig[0]-1;
	    while(i>0 && jb_dia<block_dia->Lupperneig[i-2]){
	      block_dia->Lupperneig[i] = block_dia->Lupperneig[i-2];
	      block_dia->Lupperneig[i+1] = block_dia->Lupperneig[i-1];
	      i -= 2;
	    }
	    block_dia->Lupperneig[i] = jb_dia;
	    block_dia->Lupperneig[i+1] = ioff;

	  
	  }
	  else if(jb_dia<ib_dia){

#ifdef DEBUG
	    if(block_dia->Lneig[ioff]!=jb_dia){
	  printf("error in reading list of lower neighbors in decompose!\n");
	  printf("ib_dia %d, jneig %d, jb_dia %d, ioff %d, ? %d\n",
      ib_dia, block_dia->Llowerneig[0]+1, jb_dia, ioff, block_dia->Lneig[ioff]);
	      exit(1);
	    }
	    if(block_dia->Llowerneig[0]>=99){
	      printf("increase the number of lower subdiagonal neighbors!\n");
	      exit(1);
	    }
#endif
	    
	    block_dia->Llowerneig[0]++;
	    //block_dia->Llowerneig[2*block_dia->Llowerneig[0]-1] = jb_dia;
	    //block_dia->Llowerneig[2*block_dia->Llowerneig[0]] = ioff;
	    // we must sort the array
	    i=2*block_dia->Llowerneig[0]-1;
	    while(i>0 && jb_dia<block_dia->Llowerneig[i-2]){
	      block_dia->Llowerneig[i] = block_dia->Llowerneig[i-2];
	      block_dia->Llowerneig[i+1] = block_dia->Llowerneig[i-1];
	      i -= 2;
	    }
	    block_dia->Llowerneig[i] = jb_dia;
	    block_dia->Llowerneig[i+1] = ioff;

	  }
	
	}
      }
      


/* ASSEMBLE STIFFNESS MATRIX INTO BIG BLOCKS */
      printf("Assembling preconditioner arrays\n");

/* for each diagonal block */
      for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){

/* if block is active */
	if(it_level->Block_dia[ib_dia]!=NULL){

/* to simplify typing */
	  block_dia=it_level->Block_dia[ib_dia];

/* total number of dofs for DIA */
	  nrdof_dia = block_dia->Lpos[0];

/* assemble diagonal and suitable off diagonal submatrices from
   small blocks into Dia array for big block */
/* initialize Dia and ips matrices */
	  for(i=0;i<nrdof_dia*nrdof_dia;i++) block_dia->Dia[i]=0.0;
	  for(i=0;i<nrdof_dia;i++) block_dia->Ips[i]=0;

/* loop over small blocks */
	  for(ibs=1;ibs<=block_dia->Lsmall[0];ibs++){

/* get into small block */
	    block=it_level->Block[block_dia->Lsmall[ibs]];

/* number of dofs for small block */
	    ndof1 = block->Ndof;

/* position of first dof in diagonal submatrix */
	    posloc1 = block_dia->Lpos[ibs];

/* assemble diagonal subarray */
	    for(i=0;i<ndof1;i++){
	      for(j=0;j<ndof1;j++){
		block_dia->Dia[posloc1+i+nrdof_dia*(posloc1+j)]
		  = block->Dia[i+ndof1*j];
	      }
	    }

/* if more than one small block */

	    if(block_dia->Lsmall[0]>1){

	      printf("Large blocks not working in ILU!!! Exiting\n");
	      exit(0);

/* loop over small blocks */
	      for(iaux=1;iaux<=block_dia->Lsmall[0];iaux++){

/* if different from block ibs  */
		if(iaux!=ibs){

		  ioff = it_chk_list(block_dia->Lsmall[iaux],
				     &block->Lngb[1],block->Lngb[0]);

/* if there is an off-diagonal block */
		  if(ioff){

/* number of dofs for the off-diagonal submatrix */
		    ndof2 = it_level->Block[block_dia->Lsmall[iaux]]->Ndof;

/* position for the off-diagonal submatrix in dia */
		    posloc2 = block_dia->Lpos[iaux];

/* assemble off-diagonal subarray */
		    for(i=0;i<ndof1;i++){
		      for(j=0;j<ndof2;j++){
			block_dia->Dia[posloc1+i+nrdof_dia*(posloc2+j)]
			  = block->Aux[ioff-1][i+ndof1*j];
		      }
		    }

		  } /* if there is off-diagonal block */ 
		  
		} /* if two different blocks: iaux != ibs */

	      } /* for all small blocks: iaux */

	    } /* if more small blocks */ 

	  } /* for all small blocks: ibs */

/* assemble diagonal and suitable off diagonal submatrices from
   small blocks into Aux arrays for big block */

/* initialize Aux matrices */
/* loop over big block neighbors */
	  nrneig = block_dia->Lneig[0];
	  for(ineig=0;ineig<block_dia->Lneig[0];ineig++){
	    j = block_dia->Lneig[ineig+1];
	    kaux = it_level->Block_dia[j]->Lpos[0];
	    for(i=0;i<nrdof_dia*kaux;i++) block_dia->Aux[ineig][i]=0.0;
	  }


/* loop over small blocks */
	  for(ibs=1;ibs<=block_dia->Lsmall[0];ibs++){

/* get into small block */
	    block=it_level->Block[block_dia->Lsmall[ibs]];

/* number of dofs for small block */
	    ndof1 = block->Ndof;

/* position of first dof in diagonal submatrix */
	    posloc1 = block_dia->Lpos[ibs];

/* loop over neighboring big blocks */
	    for(ineig=0;ineig<block_dia->Lneig[0];ineig++){
	      j = block_dia->Lneig[ineig+1];
	      dia_neig = it_level->Block_dia[j];

/* loop over small blocks from dia_neig */
	      for(iaux=1;iaux<=dia_neig->Lsmall[0];iaux++){
		ismneig=dia_neig->Lsmall[iaux];

/* check whether two small blocks are neighbors (non-zero entries in A) */
		ioff = it_chk_list(ismneig, &block->Lngb[1],block->Lngb[0]);

/* if there is a common off-diagonal block */
		if(ioff){

/* number of dofs for the off-diagonal submatrix */
		  ndof2 = it_level->Block[ismneig]->Ndof;

/* position for the off-diagonal submatrix in dia */
		  posloc2 = dia_neig->Lpos[iaux];

/* assemble off-diagonal subarray */
		  for(i=0;i<ndof1;i++){
		    for(j=0;j<ndof2;j++){
		      block_dia->Aux[ineig][posloc1+i+nrdof_dia*(posloc2+j)]
			= block->Aux[ioff-1][i+ndof1*j];
		    }
		  }

		} /* if there is off-diagonal block */ 
		  
	      } /* for all small blocks of a nieghbor: iaux */

	    } /* for all big neighboring blocks */ 

	  } /* for all small blocks: ibs */

	} /* if block is active */

      } /* for all blocks */


/*kbw
      //if(it_level->Monitor>ERRORS)
	{

	int nrdofbl, jaux;

	for(iblock=1;iblock<=it_level->Nrblocks;iblock++){
	  if(iblock==352){
	  block = it_level->Block[iblock];
	  nrdofbl = block->Ndof;
	  printf("Block %d, ndof %d, Posg %d\n   Neighbors:",
		 iblock,block->Ndof, block->Posg);
	  for(i=1;i<=block->Lngb[0];i++)
	    printf("  %d",block->Lngb[i]);
	  printf("\n");	  
	  printf("Dia:\n");
	  for(i=0;i<nrdofbl*nrdofbl;i++){
	    printf("%25.15lf",block->Dia[i]);
	  }
	  for(j=0;j<block->Lngb[0];j++){
	    printf("\nAux[%d]:\n",j);
	    for(i=0;i<nrdofbl*nrdofbl;i++){
	      printf("%25.15lf",block->Aux[j][i]);
	    }
	  }
	  printf("\n\n");
	  getchar();
	  }
	}

	for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){

	  //if(it_level->Block_dia[ib_dia]!=NULL){
	  if(ib_dia==352){

	    block_dia=it_level->Block_dia[ib_dia];

	    nrdofbl = block_dia->Lpos[0];
	    iaux=block_dia->Lsmall[0];
	    printf("Block_dia %d\nSmall blocks:",ib_dia);
	    for(i=1;i<=iaux;i++) printf("  %d",block_dia->Lsmall[i]);
	    printf("\n");
	    kaux=block_dia->Lneig[0];
	    printf("Neighbors:");
	    for(i=1;i<=kaux;i++) printf("  %d",block_dia->Lneig[i]);
	    printf("\n");
	    printf("Dia:\n");
	    for(i=0;i<nrdofbl*nrdofbl;i++){
	      printf("%25.15lf",block_dia->Dia[i]);
	    }
	    printf("\n");
	    nrneig = block_dia->Lneig[0];
	    for(ineig=0;ineig<block_dia->Lneig[0];ineig++){
	      j = block_dia->Lneig[ineig+1];
	      kaux = it_level->Block_dia[j]->Lpos[0];
	      printf("\nAux[%d]:\n",ineig);
	      for(i=0;i<nrdofbl*kaux;i++) 
		printf("%25.15lf",block_dia->Aux[ineig][i]);
	    }
	    printf("\n\n");
	    getchar();


	  }
	}
      }
/*kew*/

/* PERFORM INCOMPLETE LU FACTORIZATION */
      printf("Starting factorization\n");


/* for each block row (diagonal block) */
      for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){

/*kbw
	printf("block %d (%lf\%)\n",ib_dia, (double)ib_dia*100/it_level->Nrblocks_dia);
/*kew*/


/* if block is active */
	if(it_level->Block_dia[ib_dia]!=NULL){

/* to simplify typing */
	  block_dia=it_level->Block_dia[ib_dia];

/* total number of dofs for DIA */
	  nrdof_dia = block_dia->Lpos[0];

/* for all rows (blocks) already considered */
	  //for(kb_dia=1;kb_dia<ib_dia;kb_dia++){
	  //double *a_work = it_dvector(lwork,"A_work in factor_dia");
	  
	  
/* if two blocks are neighbors */
	  //ioff = it_chk_list(kb_dia,&block_dia->Lneig[1],block_dia->Lneig[0]);
	  //if(ioff){

	  int iaux_dia;
	  for(iaux_dia=1;iaux_dia<=block_dia->Llowerneig[0];iaux_dia++){

	    kb_dia = block_dia->Llowerneig[2*iaux_dia-1];
	    ioff = block_dia->Llowerneig[2*iaux_dia];

	    {


#ifdef DEBUG
	      if(block_dia->Lneig[ioff]!=kb_dia){
		printf("error in reading list of neighbors in decompose!\n");
		exit(1);
	      }
#endif

	      block_k = it_level->Block_dia[kb_dia];

/* compute product a_ik = a_ik * a_kk^-1 */
	      kaux = block_k->Lpos[0];

/*kbw
	      if(ib_dia==352){
		printf("kb_dia %d, Dia before solution:\n",kb_dia);
		for(i=0;i<kaux*kaux;i++){
		  printf("%25.15lf",block_k->Dia[i]);
		}
		printf("\n");
		printf("ib_dia %d, ioff %d, Aux before solution:\n",ib_dia,ioff);
		for(i=0;i<kaux*nrdof_dia;i++){
		  printf("%25.15lf",block_dia->Aux[ioff-1][i]);
		}
		printf("\n");
	      }
/*kew*/

              if(kaux*nrdof_dia>lwork){
		printf("Increase size of lwork to %d in factor_dia\n", 
		       kaux*nrdof_dia);
		exit(1);
	      }
	      for(i=0;i<nrdof_dia;i++){
		for(j=0;j<kaux;j++){
		  daux=0.0;
		  for(k=0;k<kaux;k++){
		    daux +=
		      block_dia->Aux[ioff-1][i+nrdof_dia*k]*
		      block_k->Dia[k+kaux*j];
		  }
		  block_dia->Aux[ioff-1][i+nrdof_dia*j]=daux;
		  //a_work[i+nrdof_dia*j] = daux;
		}
	      }
	      //for(i=0;i<nrdof_dia;i++){
	      //	for(j=0;j<kaux;j++){
	      //	  block_dia->Aux[ioff-1][i+nrdof_dia*j] = a_work[i+nrdof_dia*j];
	      //	}
	      //}
	      

/*kb!!!!!!!!!?????????????
	      dgetrs_("N",&kaux,&nrdof_dia,block_k->Dia,&kaux,block_k->Ips, 
		      block_dia->Aux[ioff-1],&kaux,&iaux);
/*kb!!!!!!!!!?????????????*/

/*kbw
	      if(ib_dia==352){
		printf("i %d (%d), k %d (%d), Aux after solution:\n",
		       ib_dia,nrdof_dia,kb_dia,kaux);
		for(i=0;i<kaux*nrdof_dia;i++){
		  printf("%25.15lf",block_dia->Aux[ioff-1][i]);
		}
		printf("\n");
	      }
/*kew*/


/* update Dia block */
	      koff = it_chk_list(ib_dia,&block_k->Lneig[1],
				 block_k->Lneig[0]);	       

/*kbw
	      if(ib_dia==352){
		printf("i %d (%d), k %d (%d), block_k->Aux for multiplication\n",
		       ib_dia,nrdof_dia,kb_dia,kaux);
		for(i=0;i<kaux*nrdof_dia;i++){
		  printf("%25.15lf",block_k->Aux[koff-1][i]);
		}
		printf("\n");
	      }
/*kew*/
	      for(i=0;i<nrdof_dia;i++){
		for(j=0;j<nrdof_dia;j++){
		  daux=0.0;
		  for(k=0;k<kaux;k++){
		    daux +=
		      block_dia->Aux[ioff-1][i+nrdof_dia*k]*
		      block_k->Aux[koff-1][k+kaux*j];
		  }
		  block_dia->Dia[i+nrdof_dia*j] -= daux;
		}
	      }

/*kbw
	      if(ib_dia==352){
		printf("ib_dia %d, Dia after update:\n",ib_dia);
		for(i=0;i<nrdof_dia*nrdof_dia;i++){
		  printf("%25.15lf",block_dia->Dia[i]);
		}
		printf("\n");
		getchar();
	      }
/*kew*/

/* loop over all neighbors of dia block, with indices greater than k */
	      nrneig = block_dia->Lneig[0];
#pragma omp parallel for default(none) private(jb_dia,dia_neig,koff,jaux,i,j,k,daux) shared(block_dia,it_level) firstprivate(kb_dia,block_k,nrdof_dia,kaux,ioff)
	      for(ineig=0;ineig<block_dia->Lneig[0];ineig++){
		jb_dia = block_dia->Lneig[ineig+1];
		if(jb_dia>kb_dia){
		  dia_neig = it_level->Block_dia[jb_dia];

/* if k-block and j-block are neighbors, i.e. Aux blocks (j,k) are non-zero */
		  koff = it_chk_list(jb_dia,&block_k->Lneig[1],
				     block_k->Lneig[0]);	       
		  if(koff){

		    jaux=dia_neig->Lpos[0];

/*kbw
  printf("i %d (%d), j %d (%d), a_ij before solution:\n",
	 ib_dia,nrdof_dia,jb_dia,jaux);
  for(i=0;i<jaux*nrdof_dia;i++){
    printf("%20.15lf",block_dia->Aux[ineig][i]);
  }
  printf("\n");
  printf("i %d (%d), k %d (%d), a_ik before solution:\n",
	 ib_dia,nrdof_dia,kb_dia,kaux);
  for(i=0;i<kaux*nrdof_dia;i++){
    printf("%20.15lf",block_dia->Aux[ioff-1][i]);
  }
  printf("\n");
  printf("k %d (%d), j %d (%d), a_kj before solution:\n",
	 kb_dia,kaux,jb_dia,jaux);
  for(i=0;i<kaux*jaux;i++){
    printf("%20.15lf",block_k->Aux[koff-1][i]);
  }
  printf("\n");
/*kew*/

/* compute a_ij = a_ij - a_ik*a_kj */

/*kb!!!!!!!!!?????????????
		    dgemm_('N','N',&nrdof_dia,&jaux,&kaux,
			   &dmone,block_dia->Aux[ioff-1],&nrdof_dia,
			   block_k->Aux[koff-1],&kaux,
			   &done,block_dia->Aux[ineig],&nrdof_dia);
/*kb!!!!!!!!!?????????????*/

		    for(i=0;i<nrdof_dia;i++){
		      for(j=0;j<jaux;j++){
			daux=0.0;
			for(k=0;k<kaux;k++){
			  daux +=
			    block_dia->Aux[ioff-1][i+nrdof_dia*k]*
			    block_k->Aux[koff-1][k+kaux*j];
			}
			block_dia->Aux[ineig][i+nrdof_dia*j] -= daux;
		      }
		    }

/*kbw
  printf("i %d, j %d, a_ij after solution:\n",ib_dia,jb_dia);
  for(i=0;i<jaux*nrdof_dia;i++){
    printf("%20.15lf",block_dia->Aux[ineig][i]);
  }
  printf("\n");
/*kew*/
/*kbw
  printf("i %d, k %d, a_ik after solution:\n",ib_dia,kb_dia);
  for(i=0;i<kaux*nrdof_dia;i++){
    printf("%20.15lf",block_dia->Aux[ioff-1][i]);
  }
  printf("\n");
  printf("k %d, j %d, a_kj after solution:\n",kb_dia,jb_dia);
  for(i=0;i<kaux*jaux;i++){
    printf("%20.15lf",block_k->Aux[koff-1][i]);
  }
  printf("\n");
  getchar();
/*kew*/

		  } /* end if non-zero entry */

		} 
	      
	      } /* end for all neighbors of dia block */
	    
	    } /* end if non-zero entry */

	    //free(a_work);	  

	  } /* end for all block already considered:  kb_dia (columns k)

/* decompose Dia */
	  if(nrdof_dia==1) block_dia->Dia[0]= 1/block_dia->Dia[0];
	  else{
/* LU decompose block->Dia */

/*kbw
	    if(ib_dia==352){
	      printf("ib_dia %d, Dia before factorization:\n",ib_dia);
	      for(i=0;i<nrdof_dia*nrdof_dia;i++){
		printf("%25.15lf",block_dia->Dia[i]);
	      }
	      printf("\n");
	    }
/*kew*/

	    double *a_work = it_dvector(lwork,"A_work in factor_dia");

	    if((nrdof_dia+1)*nrdof_dia>lwork){
	      printf("Increase size of lwork to %d in factor_dia\n", 
		     (nrdof_dia+1)*nrdof_dia);
	      exit(1);
	    }


	    it_dgetrf(block_dia->Dia,nrdof_dia,block_dia->Ips);

	    for(iaux=0;iaux<nrdof_dia;iaux++){
	      for(i=(iaux+1)*nrdof_dia;i<(iaux+2)*nrdof_dia;i++){
		a_work[i]=0.0;
	      }
	      a_work[(iaux+1)*nrdof_dia+iaux]=1.0;
	      it_dgetrs(block_dia->Dia,nrdof_dia,
			&a_work[(iaux+1)*nrdof_dia],
			&a_work[iaux*nrdof_dia],block_dia->Ips);
	    }

	    for(i=0;i<nrdof_dia*nrdof_dia;i++){
	      block_dia->Dia[i]=a_work[i];
	    }

/*kb
	    dgetrf_(&nrdof_dia,&nrdof_dia,block_dia->Dia,
		    &nrdof_dia,block_dia->Ips,&iaux);

	    dgetri_(&nrdof_dia,block_dia->Dia,
		    &nrdof_dia,block_dia->Ips,a_work,&lwork,&iaux);
/*kb*/


/*kbw
	    if(ib_dia==352){
	      
	      printf("ib_dia %d, Dia after factorization2:\n",ib_dia);
	      for(i=0;i<nrdof_dia*nrdof_dia;i++){
		printf("%25.15lf",block_dia->Dia[i]);
	      }
	      printf("\n");
	      getchar();
	    }
/*kew*/

	    free(a_work);

	  }
	      
	} /* if block is active */

      } /* for all blocks */

    

      
    } /* end if ILU preconditioning */ 

  } /* end if block size > 0 (i.e. there is Block_dia structure) */

  return(0);

}


/*---------------------------------------------------------
  it_free_blocks_dia - to free space for a block structure
---------------------------------------------------------*/
int it_free_blocks_dia(
  int Solver_id,   /* in: solver ID (used to identify the subproblem) */
  int Level_id    /* in: level ID */
  )
{

/* local variables */
  itt_levels *it_level;       /* mesh levels */
  itt_blocks_dia *block_dia;  /* to simplify typing */
  int ib_dia, ineig, nrneig;  /* counters for blocks */

/*++++++++++++++++ executable statements ++++++++++++++++*/

  
  /* get pointer to the level structure */
  it_level = &itv_solver[Solver_id].level[Level_id];
    
  for(ib_dia=1;ib_dia<=it_level->Nrblocks_dia;ib_dia++){

    block_dia = it_level->Block_dia[ib_dia];

    if( it_level->Precon==MULTI_ILU||it_level->Precon==BLOCK_ILU ){
      nrneig = block_dia->Lneig[0];
      for(ineig=0;ineig<nrneig;ineig++){
	free(block_dia->Aux[ineig]);
      }
      free(block_dia->Aux); 
    }
    
    free(block_dia->Dia);
    free(block_dia->Ips);
    
    free(block_dia->Lsmall);
    
    if(it_level->Solver==GMRES&&it_level->GMRES_type==MATRIX_FREE){
      free(block_dia->Lelem);
    }
    
    free(block_dia->Lpos);    
    
    if( it_level->Precon==MULTI_ILU||it_level->Precon==BLOCK_ILU ){
      free(block_dia->Lneig);
      free(block_dia->Llowerneig);
      free(block_dia->Lupperneig);
    }
    
    free(it_level->Block_dia[ib_dia]);
    
  } /* for all diagonal blocks */
  
  
  free(it_level->Block_dia);

  return(0);
}


/*---------------------------------------------------------
  it_free_blocks - to free space for a block structure
---------------------------------------------------------*/
int it_free_blocks(
  int Solver_id,   /* in: solver ID (used to identify the subproblem) */
  int Level_id    /* in: level ID */
  )
{

/* local variables */
  itt_levels *it_level; /* mesh levels */
  itt_blocks *block;	/* to simplify typing */

/* auxiliary variable */
  int i, iblock, bloff, nrblocks, nrdofbl, nroffbl;

/*++++++++++++++++ executable statements ++++++++++++++++*/

  /* get pointer to the level structure */
  it_level = &itv_solver[Solver_id].level[Level_id];

  /* allocate space for pointers to blocks */
  nrblocks = it_level->Nrblocks;
    
  for(iblock=1;iblock<=nrblocks;iblock++){
            
    block = it_level->Block[iblock];
     
    if(block->Lngb!=NULL){

      if(it_level->Solver!=GMRES || it_level->GMRES_type!=MATRIX_FREE){
	for(i=1;i<=block->Lngb[0];i++){
	  free(block->Aux[i-1]);
	} 
	free(block->Dia);
      }
      free(block->Aux);
      free(block->Rhs);
      free(block->Ips);
      
      free(block->Lngb);
      
    }
     
    free(it_level->Block[iblock]);
    
  }
  
  free(it_level->Block);

  return(0);
}


/*---------------------------------------------------------
it_compres - to compute the residual of the system of equations,
	v = ( b - Ax ) (used also to compute the product v = -Ax)
---------------------------------------------------------*/
void it_compres ( 
        itt_levels* It_level,/* in: pointer to current level data structure */
	int Use_rhs,	/* in: indicator whether to use RHS */
	int Ini_zero,	/* in: flag for zero initial guess */ 
	int Ndof, 	/* in: number of unknowns (components of x) */
	double* X, 	/* in: input vector (may be NULL if Ini_zero==0) */
        double* B,	/* in:  the rhs vector, if NULL take rhs */
			/*      from block data structure ( B is not taken */
	                /*      into account if Use_rhs!=1) */
	double* V 	/* out: v = b-Ax */
	)
{

/* constants */
  int ione = 1;
  double done = 1.;
  double dmone = -1.;

/* auxiliary variables */
  itt_blocks *block;  	/* to simplify typing */
  int iblock;		/* counters for blocks */
  int ioffbl;		/* counter for off diagonal blocks */
  int nrdofbl;		/* numbers of dofs in a block */
  int ndofngb;		/* number of dofs for a neighbor */

  int i,j;
  int iaux, jaux, kaux, laux, iblaux;
  double daux;

/*++++++++++++++++ executable statements ++++++++++++++++*/

/* perform matrix vector multiplication using elementary blocks */
/* loop over all small blocks */
  for(iblock=1;iblock<=It_level->Nrblocks;iblock++){

/* to simplify typing */
    block = It_level->Block[iblock];

   /* if regular (not ghost) block */
   if(block->Lngb!=NULL){

/* number of degrees of freedom for the block */
    nrdofbl = block->Ndof;

/* position of the first iblock dof in the vector v */
    kaux = block->Posg;

/* initialize part of v */
    for(i=0;i<nrdofbl;i++) V[kaux+i]=0.;

    if(!Ini_zero){
/* if initial guess not zero */

/* loop over neighbors */
      for(ioffbl=0;ioffbl<block->Lngb[0];ioffbl++){

/* global number of the subsequent neighbor */
	iaux = block->Lngb[ioffbl+1];

/* number of degrees of freedom for the neighbor */
	ndofngb = It_level->Block[iaux]->Ndof;

/* position of the first ingb dof in the global v vector */
	jaux = It_level->Block[iaux]->Posg;

/* multiplication by the off-diagonal block */
	if(nrdofbl*ndofngb>MIN_DOF_DGEMV){
	  dgemv_("N", &nrdofbl, &ndofngb, &dmone, 
		 block->Aux[ioffbl], &nrdofbl, &X[jaux], &ione, 
		 &done, &V[kaux], &ione);
	}
	else { /* if problem too small to use dgemv */
	  for(i=0;i<nrdofbl;i++){
	    for(j=0;j<ndofngb;j++){
	      V[kaux+i] -= block->Aux[ioffbl][i+j*nrdofbl]*X[jaux+j];
	    }
	  }
	}
      }

/* multiplication with diagonal block - result stored in vtemp*/
      if(nrdofbl*nrdofbl>MIN_DOF_DGEMV){
	dgemv_("N", &nrdofbl, &nrdofbl, &dmone, 
	       block->Dia, &nrdofbl, &X[kaux], &ione, 
	       &done, &V[kaux], &ione);
      }
      else { /* if problem too small to use dgemv */
	for(i=0;i<nrdofbl;i++){
	  for(j=0;j<nrdofbl;j++){
	    V[kaux+i] -= block->Dia[i+j*nrdofbl]*X[kaux+j];
	  }
	}
      }

    } /* end if not Ini_zero (not zero initial guess X) */
    
/* update using the rhs vector */
    if(Use_rhs==1){
      if(B==NULL){
	daxpy_(&nrdofbl, &done, block->Rhs, &ione, 
	       &V[kaux], &ione);
      }
      else{
	daxpy_(&nrdofbl, &done, &B[kaux], &ione, 
	       &V[kaux], &ione);
      }
    }

/*kbw
#ifdef DEBUG
    printf("compres: subsequent rows - \n");
    for(i=0;i<nrdofbl;i++) printf("%20.15lf",V[kaux+i]);
    printf("\n");
#endif
/*kew*/

   }   /* end if regular (not ghost) block */
  } /* end loop over blocks */

  return;
}


/*---------------------------------------------------------
it_blsiter - to perform one iteration of block Gauss-Seidel algorithm 
	(block Jacobi switched on by It_level->Precon==BLOCK_JACOBI)
        v_out = v_in + M^-1 * ( b - A * v_in )
---------------------------------------------------------*/
void it_blsiter(
        itt_levels* It_level,/* in: pointer to current level data structure */
	int Use_rhs,	/* in: 0 - no rhs, 1 - with rhs */
	int Ini_zero,	/* in: flag for zero initial guess */ 
	int Ndof,	/* in: number of unknowns (components of v*) */ 
	double* V,	/* in,out: vector of unknowns updated */
			/* during the loop over subdomains */
        double* B	/* in:  the rhs vector, if NULL take rhs */
			/*      from block data structure */
	)
{

/* constants */
  int ione = 1;
  double done = 1.;
  double dmone = -1.;
  
/* local variables */
  itt_blocks *block;  	/* to simplify typing */
  int iblock;		/* counter for blocks */
  int ioffbl;		/* counter for off diagonal blocks */
  int nrdofbl;		/* number of dofs in a block */
  int ndofngb;		/* number of dofs for a neighbor */
  int i_loop;		/* counter for loops for GS */
  double *vtemp;		/* temporary vector for block Jacobi*/

/* auxiliary variables */
  int i,j;
  int iaux, jaux, kaux;

/*++++++++++++++++ executable statements ++++++++++++++++*/

  if(It_level->Precon==BLOCK_JACOBI){
    vtemp = it_dvector(Ndof,"vtemp in it_bliter\n");
  }
  
/* loop over loops */
  for(i_loop=0;i_loop<It_level->Nr_prec;i_loop++){

/* loop over all blocks */
    for(iblock=1;iblock<=It_level->Nrblocks;iblock++){

/* to simplify typing */
      if(i_loop%2==0) block = It_level->Block[iblock];
      else if(i_loop%2==1) block = It_level->Block[It_level->Nrblocks+1-iblock];

/* number of degrees of freedom for the block */
      nrdofbl = block->Ndof;

/* position of the first iblock dof in the vector v */
      kaux = block->Posg;

/* initialize part of v */
      if(It_level->Precon==BLOCK_JACOBI) {
	for(i=0;i<nrdofbl;i++) vtemp[kaux+i]=0.;
      }
      else {
	for(i=0;i<nrdofbl;i++) V[kaux+i]=0.;
      }

      if(!Ini_zero||It_level->Precon==BLOCK_GS||i_loop>0){ 
/* if initial guess not zero */

/* loop over neighbors */
	for(ioffbl=0;ioffbl<block->Lngb[0];ioffbl++){

/* global number of the subsequent neighbor */
	  iaux = block->Lngb[ioffbl+1];

/* number of degrees of freedom for the neighbor */
	  ndofngb = It_level->Block[iaux]->Ndof;

/* position of the first ingb dof in the global v vector */
	  jaux = It_level->Block[iaux]->Posg;


/* solution of the local problem */
	  if(nrdofbl*ndofngb>MIN_DOF_DGEMV){
	    if(It_level->Precon==BLOCK_JACOBI) 
	      dgemv_("N", &nrdofbl, &ndofngb, &dmone, 
		     block->Aux[ioffbl], &nrdofbl, &V[jaux], &ione, 
		     &done, &vtemp[kaux], &ione);
	    else dgemv_("N", &nrdofbl, &ndofngb, &dmone, 
			block->Aux[ioffbl], &nrdofbl, &V[jaux], &ione, 
			&done, &V[kaux], &ione);
	  }
	  else { /* if problem too small to use dgemv */
	    if(It_level->Precon==BLOCK_JACOBI){
	      for(i=0;i<nrdofbl;i++){
		for(j=0;j<ndofngb;j++){
		  vtemp[kaux+i] -= block->Aux[ioffbl][i+j*nrdofbl]*V[jaux+j];
/*kbw
		  if(iblock<20){
		    printf("iblock %d, nrdofbl %d, kaux %d, ioffbl %d, iaux %d, jaux %d, i %d, j %d\n", 
			     iblock, nrdofbl, kaux, ioffbl, iaux, jaux, i, j);
		      printf("20.15lf%20.15lf%20.15lf\n", 
			     vtemp[kaux+i], block->Aux[ioffbl][i+j*nrdofbl],
			     V[jaux+j]);
		  }
/*kew*/
		}
	      }
	    }
	    else {
	      for(i=0;i<nrdofbl;i++){
		for(j=0;j<ndofngb;j++){
		  V[kaux+i] -= block->Aux[ioffbl][i+j*nrdofbl]*V[jaux+j];
		}
	      }
	    }
	  }  
	} 
      }
 
/* update using the rhs vector */
      if(Use_rhs==1){
	if(B==NULL){
	  if(It_level->Precon==BLOCK_JACOBI) {
	    daxpy_(&nrdofbl, &done, block->Rhs, &ione, 
		   &vtemp[kaux], &ione);
	  }
	  else {
	    daxpy_(&nrdofbl, &done, block->Rhs, &ione, 
		   &V[kaux], &ione);
	  }
	}
	else{
	  if(It_level->Precon==BLOCK_JACOBI) {
	    daxpy_(&nrdofbl, &done, &B[kaux], &ione, 
		   &vtemp[kaux], &ione);
	  }
	  else {
	    daxpy_(&nrdofbl, &done, &B[kaux], &ione, 
		   &V[kaux], &ione);
	  }
	}
      }
      
    } /* loop over all blocks */
    
    if(It_level->Precon==BLOCK_JACOBI){
      dcopy_(&Ndof, vtemp, &ione, V, &ione);
    }
    
  } /* loop over preconditioner iterations */
  
  if(It_level->Precon==BLOCK_JACOBI){
    free(vtemp);
  }
  
  return;
}

/*---------------------------------------------------------
it_blliter - to perform one iteration of block Gauss-Seidel algorithm 
	(block Jacobi switched on by It_level->Precon==BLOCK_JACOBI)
        v_out = v_in + M^-1 * ( b - A * v_in )
---------------------------------------------------------*/
void it_blliter(
        itt_levels* It_level,/* in: pointer to current level data structure */        
	int Use_rhs,	/* in: 0 - no rhs, 1 - with rhs */
	int Ini_zero,	/* in: flag for zero initial guess */ 
	int Ndof,	/* in: number of unknowns (components of v*) */ 
	double* V,	/* in,out: vector of unknowns updated */
			/* during the loop over subdomains */
        double* B	/* in:  the rhs vector, if NULL take rhs */
			/*      from block data structure */
	)
{

/* constants */
  int ione = 1;
  double done = 1.;
  double dmone = -1.;

/* local variables */
  itt_blocks *block;  	/* to simplify typing */
  itt_blocks_dia *block_dia;/* to simplify typing */
  int iblock, ib_dia;	/* counters for blocks */
  int ioffbl;		/* counter for off diagonal blocks */
  int nrdofbl, nrdof_dia;	/* numbers of dofs in a block */
  int ndofngb;		/* number of dofs for a neighbor */
  int i_loop;		/* counter for loops for GS */
  double vloc[Max_dof_block_dia];	/* vector for DIA dofs */
  double *vtemp;		/* temporary vector for block Jacobi*/

/* auxiliary variables */
  int i,j;
  int iaux, jaux, kaux, laux, iblaux;

/*++++++++++++++++ executable statements ++++++++++++++++*/

/*kbw
    printf("V on entrance\n");
    for(i=0;i<20;i++) printf("%20.15lf",V[i]);
    //for(i=0;i<Ndof;i++) printf("%10.6lf",V[i]);
    printf("\n");
    printf("Ini_zero %d, Use_rhs %d, Ndof %d\n", Ini_zero, Use_rhs, Ndof);
    getchar();
/*kew*/

  if(It_level->Precon==BLOCK_JACOBI){
    vtemp = it_dvector(Ndof,"vtemp in it_bliter\n");
    for(i=0;i<Ndof;i++)vtemp[i]=V[i];
  }

  for(i_loop=0;i_loop<It_level->Nr_prec;i_loop++){

/* loop over all big diagonal blocks */
    for(ib_dia=1;ib_dia<=It_level->Nrblocks_dia;ib_dia++){

      if(i_loop%2==0) iblaux = ib_dia;
      else if(i_loop%2==1) iblaux = It_level->Nrblocks_dia + 1 - ib_dia;

      if(It_level->Block_dia[iblaux]!=NULL){

/* to simplify typing */
	block_dia=It_level->Block_dia[iblaux];

/* total number of dofs for DIA */
	nrdof_dia = block_dia->Lpos[0];
	if(nrdof_dia>Max_dof_block_dia){
	  printf("Error - too much dofs in dia block, \n");
	  printf("Increase Max_dof_block_dia in itb_bliter (It_level->Block.c)\n");
	  exit(1);
	}
	
/* initialize vloc */
	for(i=0;i<nrdof_dia;i++) vloc[i]=0.;

/* first perform multiplication A*x */

/* loop over small blocks */
	for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

/* to simplify typing */
	  block = It_level->Block[block_dia->Lsmall[iblock]];

/* number of degrees of freedom for the block */
	  nrdofbl = block->Ndof;

/* position of the first iblock dof in the vector v */
	  kaux = block->Posg;

/* position in vloc */
	  laux = block_dia->Lpos[iblock];

	  if(!Ini_zero||It_level->Precon==BLOCK_GS||i_loop>0){ 
/* if initial guess not zero */

/* loop over neighbors */
	    for(ioffbl=1;ioffbl<=block->Lngb[0];ioffbl++){

/* global number of the subsequent neighbor */
	      iaux = block->Lngb[ioffbl];

/* if not found on the list of DIA blocks */
	      if(!it_chk_list(iaux,&block_dia->Lsmall[1],
			      block_dia->Lsmall[0])){

/* number of degrees of freedom for the neighbor */
		ndofngb = It_level->Block[iaux]->Ndof;

/* position of the first ingb dof in the global v vector */
		jaux = It_level->Block[iaux]->Posg;

/* multiplication by off-diagonal subarray */
		if(nrdofbl*ndofngb>MIN_DOF_DGEMV){
		  dgemv_("N", &nrdofbl, &ndofngb, &dmone, 
			 block->Aux[ioffbl-1], &nrdofbl, &V[jaux], &ione, 
			 &done, &vloc[laux], &ione);
		}
		else { /* if problem too small to use dgemv */
		  for(i=0;i<nrdofbl;i++){
		    for(j=0;j<ndofngb;j++){
		      vloc[laux+i] -=block->Aux[ioffbl-1][i+j*nrdofbl]*V[jaux+j];
/*kbw
		      if(iblaux<20){
		      printf("%6d%6d%6d%20.15lf%20.15lf%20.15lf\n", 
			     iblaux, block_dia->Lsmall[iblock], ioffbl,
			     vloc[laux+i], block->Aux[ioffbl-1][i+j*nrdofbl],
			     V[jaux+j]);
		      }
/*kew*/
		    }
		  }
		}
/*kbw
    getchar();
/*kew*/
		
	      } /* if not found on the list of DIA blocks */
	    } /* loop over neighbors */
	  } /* if initial guess not zero */
	  
/*kbw
	  if(iblaux<20){
    printf("V after multiplication\n");
    for(i=0;i<nrdof_dia;i++) printf("%10.6lf",vloc[i]);
    printf("\n");
    //getchar();
	  }
/*kew*/

/* update using the rhs vector */
	  if(Use_rhs==1){
	    if(B==NULL){
	      daxpy_(&nrdofbl, &done, block->Rhs, &ione, 
		     &vloc[laux], &ione);
	    }
	    else{
	      daxpy_(&nrdofbl, &done, &B[kaux], &ione, 
		     &vloc[laux], &ione);
	    }


/*kbw
	  if(iblaux<20){
    printf("V after RHS subtraction\n");
    for(i=0;i<nrdof_dia;i++) printf("%10.6lf",vloc[i]);
    printf("\n");
    //getchar();
	  }
/*kew*/

	  } /* if updated using the rhs vector */

	} /* loop over small blocks */


/* resolve subdomain problem - solve for diagonal submatrix */
/*   v = block->Dia^(-1) * vloc */
       
	if(nrdof_dia==1){
	  vloc[0]*=block_dia->Dia[0];
	}
	else{
	  dgetrs_("N",&nrdof_dia,&ione,block_dia->Dia,&nrdof_dia,
		  block_dia->Ips,vloc,&nrdof_dia,&iaux);
	}

/* rewrite back to global vector */
/* loop over small blocks */
	for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

/* to simplify typing */
	  block = It_level->Block[block_dia->Lsmall[iblock]];

/* number of degrees of freedom for the block */
	  nrdofbl = block->Ndof;

/* position of the first iblock dof in the vector v */
	  kaux = block->Posg;

/* position in vloc */
	  laux = block_dia->Lpos[iblock];
	  
	  if(It_level->Precon==BLOCK_JACOBI) {
/* assume each small block is updated only once */

	    for(i=0;i<nrdofbl;i++) vtemp[kaux+i] = vloc[laux+i];
	    
	  }
	  else {
	    for(i=0;i<nrdofbl;i++) V[kaux+i] = vloc[laux+i];
	  }

	} /* loop over small blocks */
	
      } /* if block_dia is active */
    } /* loop over big diagonal blocks */
    
    if(It_level->Precon==BLOCK_JACOBI){
      dcopy_(&Ndof, vtemp, &ione, V, &ione);
    }
    
  } /* loop over preconditioner iterations */


  if(It_level->Precon==BLOCK_JACOBI){
    free(vtemp);
  }

/*kbw
    printf("V on leaving\n");
    for(i=0;i<20;i++) printf("%20.15lf",V[i]);
    //for(i=0;i<Ndof;i++) printf("%10.6lf",V[i]);
    printf("\n");
    getchar();
/*kew*/

  return;
}

/*---------------------------------------------------------
it_rhsub - to perform forward reduction and back-substitution for ILU
           preconditioning: v_out = M^-1 * b
---------------------------------------------------------*/
void it_rhsub(
        itt_levels* It_level,/* in: pointer to current level data structure */
	int Ndof,	/* in: number of unknowns (components of v*) */ 
	double* V,	/* out: vector of unknowns updated */
			/*      during the loop over subdomains */
        double* B	/* in:  the rhs vector, if NULL take rhs */
			/*      from block data structure */
	)
{

/* constants */
  int ione = 1;
  double done = 1.;
  double dmone = -1.;


/* local variables */
  itt_blocks *block, *blneig;  	/* to simplify typing */
  itt_blocks_dia *block_dia, *blneig_dia;/* to simplify typing */
  int iblock, ioffbl, ib_dia, jb_dia, ibl_sm, jbl_sm; /* counters for blocks */
  int nrdof1, nrdof2, nrdof_dia;	/* numbers of dofs in a block */
  int posloc1, posloc2, posglob1, posglob2; /* positions in dofs vectors */
  double vloc[Max_dof_block_dia];	/* vector for DIA dofs */
  double vloc1[Max_dof_block_dia];	/* vector for DIA dofs */
 
/* auxiliary variables */
  int i,j, jneig, iaux, ioff;

/*++++++++++++++++ executable statements ++++++++++++++++*/


/* FORWARD REDUCTION */

/* loop over all big diagonal blocks */
  for(ib_dia=1;ib_dia<=It_level->Nrblocks_dia;ib_dia++){

/* to simplify typing */
    block_dia=It_level->Block_dia[ib_dia];

/* total number of dofs for DIA */
    nrdof_dia = block_dia->Lpos[0];
    if(nrdof_dia>Max_dof_block_dia){
      printf("Error - too much dofs in dia block, \n");
      printf("Increase Max_dof_block_dia in itb_bliter (It_level->Block.c)\n");
      exit(1);
    }
	
/* initialize vloc */
    for(i=0;i<nrdof_dia;i++) vloc[i]=0.;

/* substitute value of rhs to vloc */

/* loop over small blocks */
    for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

      ibl_sm = block_dia->Lsmall[iblock];

/* to simplify typing */
      block = It_level->Block[ibl_sm];

/* number of degrees of freedom for the block */
      nrdof1 = block->Ndof;

/* position of the first iblock dof in the vector v */
      posglob1 = block->Posg;

/* position in vloc */
      posloc1 = block_dia->Lpos[iblock];

/* update using the rhs vector */
      if(B==NULL){
	daxpy_(&nrdof1, &done, block->Rhs, &ione, 
	       &vloc[posloc1], &ione);
      }
      else{
	daxpy_(&nrdof1, &done, &B[posglob1], &ione, 
	       &vloc[posloc1], &ione);
      }

/* for all neighboring factorized (lower subdiagonal) entries of A */
      for(jneig=1;jneig<=block_dia->Llowerneig[0];jneig++){

	jb_dia = block_dia->Llowerneig[2*jneig-1];
	ioff = block_dia->Llowerneig[2*jneig];

#ifdef DEBUG
	  if(block_dia->Lneig[ioff]!=jb_dia){
	    printf("error in reading list of lower neighbors in rhsub!\n");
	    printf("ib_dia %d, iblock %d, jneig %d, jb_dia %d, ioff %d, ? %d\n",
		   ib_dia, iblock, jneig, jb_dia, ioff, block_dia->Lneig[ioff]);
	    exit(1);
	  }
#endif

	  blneig_dia = It_level->Block_dia[jb_dia];


/* loop over neighbor's small blocks */
	  for(ioffbl=1;ioffbl<=blneig_dia->Lsmall[0];ioffbl++){

/* global number of the subsequent neighbor */
	    jbl_sm = blneig_dia->Lsmall[ioffbl];
	    blneig = It_level->Block[jbl_sm];

/* number of degrees of freedom for the neighbor */
	      nrdof2 = blneig->Ndof;

/* position of the first ingb dof in the global v vector */
	      posglob2 = blneig->Posg;

/* position of the first ingb dof in the Aux array  */
	      posloc2 = blneig_dia->Lpos[ioffbl];

/* multiplication by off-diagonal subarray 
		if(nrdofbl*ndofngb>MIN_DOF_DGEMV){
		  dgemv_("N", &nrdofbl, &ndofngb, &dmone, 
			 block->Aux[ioffbl-1], &nrdofbl, &V[jaux], &ione, 
			 &done, &vloc[laux], &ione);
		}
		else {
*/		
	      for(i=0;i<nrdof1;i++){
		for(j=0;j<nrdof2;j++){
		  iaux = posloc1+i+(posloc2+j)*nrdof_dia;
		  vloc[posloc1+i] -= 
		    block_dia->Aux[ioff-1][iaux]*V[posglob2+j];
		}
	      }

/*		} */
		
/*kbw
  printf("ib_dia %d, jb_dia %d, Aux before solution:\n",ib_dia,jb_dia);
  for(i=0;i<blneig_dia->Lpos[0]*nrdof_dia;i++){
    printf("%20.15lf",block_dia->Aux[ioff-1][i]);
  }
  printf("\n");
printf("\nblock_dia %d small block %d after updating with small block %d\n",
       ib_dia,ibl_sm,jbl_sm);
for(i=0;i<nrdof_dia;i++) printf("%10.6lf",vloc[i]);
getchar();
/*kew*/

	  } /* loop over small blocks of a neighboring big blocks */
      } /* loop over neighboring big blocks */
    } /* loop over small blocks */
   

/*kbw
  printf("kb_dia %d, Dia before solution:\n",kb_dia);
  for(i=0;i<kaux*kaux;i++){
    printf("%20.15lf",block_k->Dia[i]);
  }
  printf("\n");
  printf("ib_dia %d, ioff %d, Aux before solution:\n",ib_dia,ioff);
  for(i=0;i<kaux*nrdof_dia;i++){
    printf("%20.15lf",block_dia->Aux[ioff-1][i]);
  }
  printf("\n");
/*kew*/

	      

/* rewrite back to global vector */
/* loop over small blocks */
    for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

/* to simplify typing */
      block = It_level->Block[block_dia->Lsmall[iblock]];

/* number of degrees of freedom for the block */
      nrdof1 = block->Ndof;

/* position of the first iblock dof in the vector v */
      posglob1 = block->Posg;

/* position in vloc */
      posloc1 = block_dia->Lpos[iblock];
	  
      for(i=0;i<nrdof1;i++) V[posglob1+i] = vloc[posloc1+i];

    } /* loop over small blocks */
	
  } /* loop over big diagonal blocks */
    
/* BACK SUBSTITUTION */
/* loop over all big diagonal blocks */
  for(ib_dia=It_level->Nrblocks_dia;ib_dia>=1;ib_dia--){

/* to simplify typing */
    block_dia=It_level->Block_dia[ib_dia];

/* total number of dofs for DIA */
    nrdof_dia = block_dia->Lpos[0];
	

/* loop over small blocks */
    for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

      ibl_sm = block_dia->Lsmall[iblock];

/* to simplify typing */
      block = It_level->Block[ibl_sm];

/* number of degrees of freedom for the block */
      nrdof1 = block->Ndof;

/* position of the first iblock dof in the vector v */
      posglob1 = block->Posg;

/* position in vloc */
      posloc1 = block_dia->Lpos[iblock];

/* initialize vloc */
      for(i=0;i<nrdof1;i++) vloc[posloc1+i] = V[posglob1+i] ;

/* for all neighboring factorized (upper subdiagonal) entries of A */
      for(jneig=1;jneig<=block_dia->Lupperneig[0];jneig++){

	jb_dia = block_dia->Lupperneig[2*jneig-1];
	ioff = block_dia->Lupperneig[2*jneig];

#ifdef DEBUG
	if(block_dia->Lneig[ioff]!=jb_dia){
	  printf("error in reading list of upper neighbors in rhsub!\n");
	  printf("ib_dia %d, iblock %d, jneig %d, jb_dia %d, ioff %d, ? %d\n",
		 ib_dia, iblock, jneig, jb_dia, ioff, block_dia->Lneig[ioff]);
	  exit(1);
	}
#endif

	  blneig_dia = It_level->Block_dia[jb_dia];

/* loop over neighbor's small blocks */
	  for(ioffbl=1;ioffbl<=blneig_dia->Lsmall[0];ioffbl++){

/* global number of the subsequent neighbor */
	    jbl_sm = blneig_dia->Lsmall[ioffbl];
	    blneig = It_level->Block[jbl_sm];

/* number of degrees of freedom for the neighbor */
	      nrdof2 = blneig->Ndof;

/* position of the first ingb dof in the global v vector */
	      posglob2 = blneig->Posg;

/* position of the first ingb dof in the Aux array  */
	      posloc2 = blneig_dia->Lpos[ioffbl];

/* multiplication by off-diagonal subarray 
		if(nrdofbl*ndofngb>MIN_DOF_DGEMV){
		  dgemv_("N", &nrdofbl, &ndofngb, &dmone, 
			 block->Aux[ioffbl-1], &nrdofbl, &V[jaux], &ione, 
			 &done, &vloc[laux], &ione);
		}
		else {
*/		
	      for(i=0;i<nrdof1;i++){
		for(j=0;j<nrdof2;j++){
		  iaux = posloc1+i+(posloc2+j)*nrdof_dia;
		  vloc[posloc1+i] -= 
		    block_dia->Aux[ioff-1][iaux]*V[posglob2+j];
		}
	      }

/*		} */
		
/*kbw
  printf("ib_dia %d, jb_dia %d, Aux after solution:\n",ib_dia,jb_dia);
  for(i=0;i<blneig_dia->Lpos[0]*nrdof_dia;i++){
    printf("%20.15lf",block_dia->Aux[ioff-1][i]);
  }
  printf("\n");
printf("\nblock_dia %d small block %d after updating with small block %d\n",
       ib_dia,ibl_sm,jbl_sm);
for(i=0;i<nrdof_dia;i++) printf("%10.6lf",vloc[i]);
getchar();
/*kew*/

	  } /* loop over small blocks of a neighboring big blocks */
      } /* loop over neighboring big blocks */
    } /* loop over small blocks */
   
/* resolve subdomain problem - solve for diagonal submatrix */
/*   v = block->Dia^(-1) * vloc */
       
    if(nrdof_dia==1){
      vloc1[0]=block_dia->Dia[0]*vloc[0];
    }
    else{
	      for(i=0;i<nrdof_dia;i++){
		  vloc1[i]=0.0;
		for(j=0;j<nrdof_dia;j++){
		    vloc1[i] +=
		      block_dia->Dia[i+nrdof_dia*j]*vloc[j];
		  }
	      }

    }

/*kbw
printf("\nblock_dia %d small block %d before solving with Dia\n",
       ib_dia,ibl_sm);
for(i=0;i<nrdof_dia;i++) printf("%10.6lf",vloc[i]);
  printf("ib_dia %d, Dia before solution:\n",ib_dia);
  for(i=0;i<nrdof_dia*nrdof_dia;i++){
    printf("%20.15lf",block_dia->Dia[i]);
  }
  printf("\n");
printf("\nblock_dia %d small block %d after solving with Dia\n",
       ib_dia,ibl_sm);
for(i=0;i<nrdof_dia;i++) printf("%10.6lf",vloc1[i]);
getchar();
/*kew*/


/* rewrite back to global vector */
/* loop over small blocks */
    for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

/* to simplify typing */
      block = It_level->Block[block_dia->Lsmall[iblock]];

/* number of degrees of freedom for the block */
      nrdof1 = block->Ndof;

/* position of the first iblock dof in the vector v */
      posglob1 = block->Posg;

/* position in vloc */
      posloc1 = block_dia->Lpos[iblock];
	  
      for(i=0;i<nrdof1;i++) V[posglob1+i] = vloc1[posloc1+i];

    } /* loop over small blocks */
	
  } /* loop over big diagonal blocks */
    
  return;
}


/*---------------------------------------------------------
it_mfaiter - to perform matrix vector multiplication (possibly
	in matrix-free manner) and additive Schwarz approximate solve
---------------------------------------------------------*/
void it_mfaiter(
        itt_levels* It_level, /* in: pointer to current level data structure */
	int Use_rhs,	/* in: 0 - no rhs, 1 - with rhs */
	int Ini_zero,	/* in: flag for zero initial guess */ 
	int Ndof,	/* in: number of unknowns (components of v*) */ 
	double* V,	/* in,out: vector of unknowns updated */
			/* during the loop over subdomains */
        double* B	/* in:  the rhs vector, if NULL take rhs */
			/*      from block data structure */
	)
{

/* constants */
int ione = 1;
double done = 1.;
double dmone = -1.;

/* local variables */
itt_blocks *block;  	/* to simplify typing */
itt_blocks_dia *block_dia;/* to simplify typing */
int iblock, ib_dia;	/* counters for blocks */
int ioffbl;		/* counter for off diagonal blocks */
int nrdofbl, nrdof_dia;	/* numbers of dofs in a block */
int ndofngb;		/* number of dofs for a neighbor */
int i_loop;		/* counter for loops for GS */
double vloc[Max_dof_block_dia];	/* vector for DIA dofs */
double inveps;	/* inverse of small parameter for finite differences */

/* auxiliary variables */
int i,j;
int iaux, jaux, kaux, laux, iblaux, iel;
double daux, *vtemp, *vtemp1;

/*++++++++++++++++ executable statements ++++++++++++++++*/

/*kb!!!!!! inveps = 1.0/It_level->Jaceps; */

vtemp = it_dvector(Ndof,"vtemp in itb_dditer\n");
vtemp1 = it_dvector(Ndof,"vtemp in itb_dditer\n");

/* loop over loops */
for(i_loop=0;i_loop<It_level->Nr_prec;i_loop++){

/* first perform matrix vector multiplication using elementary blocks */
/* make global pointer It_level->u_MF to point to current vector v 
It_level->u_MF=V; */

/* for each element active in a given mesh*/
  //for(iel=1;iel<=It_level->Mesh_elems[0];iel++){
 
/* compute element load vector for matrix free multiplication */
/* and assemble to block arrays Aux */
/*  pbr_assemble_matrix_free(iel); */

  //}

/* loop over all blocks */
for(iblock=1;iblock<=It_level->Nrblocks;iblock++){

/* to simplify typing */
  if(i_loop%2==0) block = It_level->Block[iblock];
  else if(i_loop%2==1) block = It_level->Block[It_level->Nrblocks+1-iblock];

/* number of degrees of freedom for the block */
  nrdofbl = block->Ndof;

/* position of the first iblock dof in the vector v */
  kaux = block->Posg;

/* product Ax in a matrix free manner */
  if(B==NULL){
    for(i=0;i<nrdofbl;i++){
      vtemp[kaux+i] = inveps*(block->Aux[0][i]-block->Rhs[i]);
    }
  }
/* for matrix free formulation B may store modified RHS */
  else{
    for(i=0;i<nrdofbl;i++){
      vtemp[kaux+i] = inveps*(block->Aux[0][i]-block->Rhs[i]);
    }
  }

/* update using the rhs vector */
  if(Use_rhs==1){
    if(B==NULL){
      daxpy_(&nrdofbl, &done, block->Rhs, &ione, 
        &vtemp[kaux], &ione);
    }
/* for matrix free formulation B may store modified RHS */
    else{
      daxpy_(&nrdofbl, &done, block->Rhs, &ione, 
        &vtemp[kaux], &ione);
    }
  }

} /* loop over all blocks */

/* clear Aux arrays for new iterations */
for(iblock=1;iblock<=It_level->Nrblocks;iblock++){
  for(iaux=0;iaux<It_level->Block[iblock]->Ndof;iaux++){
    It_level->Block[iblock]->Aux[0][iaux]=0;
  } 
}

/* create vector for initial state before solving subdomain problems */
dcopy_(&Ndof, vtemp, &ione, vtemp1, &ione);

/* perform preconditioning step using big blocks */

/* loop over all big diagonal blocks */
for(ib_dia=1;ib_dia<=It_level->Nrblocks_dia;ib_dia++){

  if(i_loop%2==0) iblaux = ib_dia;
  else if(i_loop%2==1) iblaux = It_level->Nrblocks_dia + 1 - ib_dia;

  if(It_level->Block_dia[iblaux]!=NULL){

/* to simplify typing */
    block_dia=It_level->Block_dia[iblaux];

/* total number of dofs for DIA */
    nrdof_dia = block_dia->Lpos[0];
    if(nrdof_dia>Max_dof_block_dia){
      printf("Error - too much dofs in dia block, \n");
      printf("Increase Max_dof_block_dia in itb_bliter (It_level->Block.c)\n");
      exit(1);
    }

/* loop over small blocks */
    for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

/* to simplify typing */
      block = It_level->Block[block_dia->Lsmall[iblock]];

/* number of degrees of freedom for the block */
      nrdofbl = block->Ndof;

/* position of the first iblock dof in the vector v */
      kaux = block->Posg;

/* position in vloc */
      laux = block_dia->Lpos[iblock];

/* rewrite vtemp to vloc */
      for(i=0;i<nrdofbl;i++) vloc[laux+i] = vtemp1[kaux+i];

    } /* loop over small blocks */

/* resolve subdomain problem - solve for diagonal submatrix */
/*   vloc = block->Dia^(-1) * vloc */
       
    if(nrdof_dia==1){
      vloc[0]*=block_dia->Dia[0];
    }
    else{
      dgetrs_("N",&nrdof_dia,&ione,block_dia->Dia,&nrdof_dia,
	      block_dia->Ips,vloc,&nrdof_dia,&iaux);
    }

/* rewrite back to global vector */
/* loop over small blocks */
    for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

/* to simplify typing */
      block = It_level->Block[block_dia->Lsmall[iblock]];

/* number of degrees of freedom for the block */
      nrdofbl = block->Ndof;

/* position of the first iblock dof in the vector v */
      kaux = block->Posg;

/* position in vloc */
      laux = block_dia->Lpos[iblock];

/* we use the newest values in vloc to update vtemp */
      for(i=0;i<nrdofbl;i++) vtemp[kaux+i] = vloc[laux+i];

    } /* loop over small blocks */

  } /* if block_dia is active */
} /* loop over big diagonal blocks */

/* add v to get -A_ii^-1 * SUM_j ( A_ij * x_j - b_i ) (j!=i) */
daux = 1;
daxpy_(&Ndof, &daux, vtemp, &ione, V, &ione);

} /* loop over standard iterations in preconditioning */

free(vtemp);
free(vtemp1);

return;
}

/*---------------------------------------------------------
it_mfmiter - to perform one iteration of block Gauss-Seidel
	algorithm (block Jacobi switched on by It_level->Precon==BLOCK_JACOBI)
	within Newton-Krylov-Schwarz algorithm
---------------------------------------------------------*/
void it_mfmiter(
        itt_levels* It_level,/* in: pointer to current level data structure */        
	int Use_rhs,	/* in: 0 - no rhs, 1 - with rhs */
	int Ini_zero,	/* in: flag for zero initial guess */ 
	int Ndof,	/* in: number of unknowns (components of v*) */ 
	double* V,	/* in,out: vector of unknowns updated */
			/* during the loop over subdomains */
        double* B	/* in:  the rhs vector, if NULL take rhs */
			/*      from block data structure */
	)
{

/* constants */
int ione = 1;
double done = 1.;

/* local variables */
itt_blocks *block;  	/* to simplify typing */
itt_blocks_dia *block_dia;/* to simplify typing */
int iblock, ib_dia;	/* counters for blocks */
int nrdofbl, nrdof_dia;	/* numbers of dofs in a block */
int i_loop;		/* counter for loops for GS */
double vloc[Max_dof_block_dia];	/* vector for DIA dofs */
double inveps;	/* inverse of small parameter for finite differences */

/* auxiliary variables */
int i,iel;
int iaux, kaux, laux, iblaux;

/*++++++++++++++++ executable statements ++++++++++++++++*/

/* inveps = 1.0/it_Jaceps; */

/* make global pointer it_u_MF to point to current vector v 
it_u_MF=V; */

for(i_loop=0;i_loop<It_level->Nr_prec;i_loop++){

/* loop over all big diagonal blocks */
for(ib_dia=1;ib_dia<=It_level->Nrblocks_dia;ib_dia++){

  if(i_loop%2==0) iblaux = ib_dia;
  else if(i_loop%2==1) iblaux = It_level->Nrblocks_dia + 1 - ib_dia;

  if(It_level->Block_dia[iblaux]!=NULL){

/* to simplify typing */
    block_dia=It_level->Block_dia[iblaux];

/* total number of dofs for DIA */
    nrdof_dia = block_dia->Lpos[0];
    if(nrdof_dia>Max_dof_block_dia){
      printf("Error - too much dofs in dia block, \n");
      printf("Increase Max_dof_block_dia in itb_bliter (It_level->Block.c)\n");
      exit(1);
    }

/* first perform multiplication A*x */

/* loop over all subdomain elements */
    for(iel=1;iel<=block_dia->Lelem[0];iel++){
/*      pbr_assemble_matrix_free(block_dia->Lelem[iel]); */
    }


/* loop over small blocks */
    for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

/* to simplify typing */
      block = It_level->Block[block_dia->Lsmall[iblock]];

/* number of degrees of freedom for the block */
      nrdofbl = block->Ndof;

/* position of the first iblock dof in the vector v */
      kaux = block->Posg;

/* position in vloc */
      laux = block_dia->Lpos[iblock];

/* product Ax in a matrix free manner */
      if(B==NULL){
        for(i=0;i<nrdofbl;i++){
          vloc[laux+i] = inveps*(block->Aux[0][i]-block->Rhs[i]);
        } 
      }
      else{
/* for matrix free formulation B may store modified RHS */
        for(i=0;i<nrdofbl;i++){
          vloc[laux+i] = inveps*(block->Aux[0][i]-block->Rhs[i]);
        } 
      }

/* update using the rhs vector */
      if(Use_rhs==1){
        if(B==NULL){
          daxpy_(&nrdofbl, &done, block->Rhs, &ione, 
	        &vloc[laux], &ione);
        }
        else{
/* for matrix free formulation B may store modified RHS */
          daxpy_(&nrdofbl, &done, block->Rhs, &ione, 
	        &vloc[laux], &ione);
        }

      } /* if updated using the rhs vector */

    } /* loop over small blocks */

/* resolve subdomain problem - solve for diagonal submatrix */
/*   v = block->Dia^(-1) * vloc */
       
    if(nrdof_dia==1){
      vloc[0]*=block_dia->Dia[0];
    }
    else{
      dgetrs_("N",&nrdof_dia,&ione,block_dia->Dia,&nrdof_dia,
	      block_dia->Ips,vloc,&nrdof_dia,&iaux);
    }

/* rewrite back to global vector */
/* loop over small blocks */
    for(iblock=1;iblock<=block_dia->Lsmall[0];iblock++){

/* to simplify typing */
      block = It_level->Block[block_dia->Lsmall[iblock]];

/* number of degrees of freedom for the block */
      nrdofbl = block->Ndof;

/* position of the first iblock dof in the vector v */
      kaux = block->Posg;

/* position in vloc */
      laux = block_dia->Lpos[iblock];

      for(i=0;i<nrdofbl;i++) V[kaux+i] += vloc[laux+i];

    } /* loop over small blocks */

/* clean Aux arrays */
    for(iel=1;iel<=block_dia->Lelem[0];iel++){
/*      pbr_assemble_matrix_free(-block_dia->Lelem[iel]); */
    }

  } /* if block_dia is active */
} /* loop over big diagonal blocks */

} /* loop over preconditioner iterations */

return;
}

